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ABSTRACT:   Recently  the  solution  of  some  nonlinear  free 
boundary  problems  has  been  effected  niimerically  by 
incorporating  a   change  of  independent  variables  through 
a  conformal  transformation,  thus  simplifying  the  domain  of 
the  solution  (cf.  E.  Bloch,  [1]).   One  then  solves  the 
transformed  differential  equation  in  the  simpler  domain, 
simultaneously  determining  the  transformation  itself,  and 
the  result  is  regarded  as  a  parametrized  form  of  the  desired 
solution.   There  is  some  question  as  to  what  is  the  best 
way  to  handle  the  boundary  conditions  in  such  an  approach . 
The  aforementioned  report  employed  a  steepest  descent 
procedure  to  produce  difference  equations  at  the  boundary 
which  uniquely  determined  the  solution,  but  led  to  a  rather 
slow  iterative  scheme. 

The  present  paper  discusses  the  results  of  using 
the  reflection  property  of  solutions  of  elliptic  partial 
differential  equations  [cf.  8]  to  determine  the  boundary 
conditions  to  the  transformed  difference  equations.  The 
basic  idea  is  analyzed  theoretically,  and  we  demonstrate 
its  applicability  to  a  special  class  of  two-dimensional 
problems.   The  procedure  is  then  applied  to  a  simple 
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plasma  containment  problem  and  to  the  two-dimensional  and 
the  three-dimensional  axially  symmetric  vena  contracta 
models.   The  results  are  compared  with  those  of  Bloch, 
and  it  is  seen  that  the  reflection  scheme  requires  about 
one-tenth  as  many  iterations  for  convergence  as  the  method 
of  steepest  descent.   New  computations  for  the   contraction 
coefficient  are  presented. 
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Part  1.   INTRODUCTION 

The  advent  of  the  modern  computer  has  stimulated 
the  finite  difference  technique  of  finding  approximations 
to  solutions  of  partial  differential  equations.   Nonetheless 
the  enormous  number  of  calculations  required  by  this 
technique  strains  the  capabilities  of  even  the  biggest 
machines,  necessitating  studies  seeking  more  efficient  ways 
of  formulating  and  solving  the  equations .   In  the  case  of 
elliptic  boundary  value  problems  one  is  usually  confronted 
with  the  task  of  solving  a  large  system  of  algebraic  equa- 
tions expressing  the  finite  difference  analogues  of  the 
partial  differential  equation  and  the  boundary  conditions. 
The  solution  of  this  system  is  achieved  by  some  iterative 
procedure.   Thus  the  effectiveness  and  feasibility  of 
a  finite  difference  scheme  are  indicated, by  first,  the 
rate  at  which  the  iterates  of  the  trial  solutions  of  the 
difference  equations  converge  to  the  actual  solution, 
and  second,  the  degree  to  which  this  solution  of  the 
difference  equations  approximates  the  solution  of  the 
differential  equation.   The  latter  is  a  measure  of  the 
accuracy  of  the  finite  difference  approximation. 

The  accuracy  and  rate  of  convergence  of  many 
linear  systems  can  be  estimated  theoretically  [2,4,6,10]. 
The  appearance  of  nonlinearities ,  however,  makes  these 
estimates  prohibitively  difficult.   In  fact,  in  the  case 
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of  nonlinear  boundary  conditions ,  particularly  those 
involving  free  boundaries,  it  may  not  be  clear  how  to 
formulate  the  finite  difference  approximation.   In  this 
paper  we  describe  a  method  for  handling  some  of  these 
nonlinear  examples .   Applying  the  technique  to  a  free 
bounaary  problem,  we  are  also  able  to  establish  the 
accuracy  and  convergence  of  our  scheme  by  experimenting 
with  tne  number  of  mesh  points. 

In  the  past,  free  boundary  problems  have  been 
attacked  by  iterative  procedures  which  guess  the  location 
of  the  boundary,  solve  the  rest  of  the  equations,  and  then 
improve  the  guess  [10] .   A  more  accurate  and  faster  scheme 
was  recently  reported  by  Bloch  [1].   He  incorporates  a  change 
of  independent  variables  via  a  conformal  transformation, 
simplifying  the  domain  of  the  solution.   Then  he  solves  the 
transformea  differential  equation  in  the  new  (known)  domain, 
simultaneously  solving  for  the  transformation  itself.   The 
result  is  regarded  as  a  parametrized   form  of  the  desired 
solution.   The  boundary  conditions  for  the  transformation 
are  derived  from  a  steepest  descent  argioment,  supplemented 
by  equations  serving  to  establish  the  free  boundary 
constraint,  which  in  his  case  expresses  continuity  of 
pressure  across  a  fluid-air  interface. 

In  this  paper  we  propose  an  improvement  of  this 
procedure  whereby  the  boundary  conditions  of  the   transforma- 
tion are  obtained  by  analytic  continuation,  resulting  in 
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reflection  rules  (cf.  [8])  replacing  the  steepest  descent 
and  constant  pressure  equations  of  the  above  method.   The 
scheme  is  just  as  accurate  as  Bloch's,  but  the  iterations 
converge  almost  ten  times  faster.   Furthermore,  it  is 
directly  applicable  to  any  problem  in  which  the  solution 
can  be  continued  across  the  boundary  in  some  prescribed 
manner  as  described  by  H.  Lewy  in  [8] ,  for  example.   This 
report  purports  to  demonstrate  that  the  reflection  technique 
is  superior  to  the  method  of  steepest  descent  in  handling 
these  problems . 

In  Part  2  of  this  paper  we  describe  the  theory  behind 
the  reflection  scheme.   It  is  introduced  in  the  context  of  a 
general  free  surface  problem  of  fluid  mechanics,  but  the 
actual  procedure  for  obtaining  the  conformal  map  is  presented 
without  such  motivation.   Proofs  of  certain  results  are 
given  when  available,  and  evidence  of  other  expected 
consequences  is  offered  so  that  a  heuristic  understanding 
of  the  scheme's  implications  and  its  limitations  is  obtained. 
In  fact,  the  whole  theory  can  be  neatly  capsulized,  and 
this  is  done  in  a  summary  at  the  end  of  the  section. 

The  application  of  the  procedure  in  solving  the 
vena  contracta  problem  is  given  in  Part  3.   The  physical 
situation  is  the  following:   one  has  an  incompressible 
inviscid  fluid  under  high  pressure  confined  behind  a  plane 
wall,  and  the  wall  has  an  aperture  through  which  the  fluid 
escapes  as  a  jet.   In  the  two-dimensional  model  the 
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aperture  is  an  infinitely  long  slit;  in  the  three-dimensional 
proDlem  with  axial  symmetry,  the  aperture  is  a  circular  hole. 
The  boundary  of  the  jet  is  of  course  unknown  beforehand, 
so  we  have  to  deal  with  a  free  boundary  problem. 

To  determine  the  vena  contracta  one  must  modify 
tne  general  procedure  of  Part  2  somewhat,  so  additional  study 
of  tne  theoretical  aspects  is  necessary.   Furthermore,  this 
problem  is  full  of  complications  not  related  to  reflection, 
and  each  of  these  must  be  dealt  with.   In  particular,  the 
behavior  of  the  flow  near  the  edge  of  the  aperture  is  quite 
singular,  and  we  develop  some  special  techniques  for  handling 
this  point.   Knowing  the  location  of  this  edge,  one  can 
evaluate  the  contraction  coefficient,  which  is  the  ratio 
of  the  area  of  the  jet  at  infinity  to  the  area  of  the 
aperture.   Tnus  we  are  able  to  present  a  new,  accurate  method 
for  computing  this  number.   The  details  of  all  these 
considerations  are  presented  in  Part  3,  which  concludes 
with  a  description  of  a  Fortran  program  (listed  in  Appendix  II) 
for  the  solution  of  the  problem. 

Part  4  describes  the  results  of  running  this  program 
on  New  York  University's  CDC  6600  computer.   The  accuracy 
and  efficiency  are  reported,  and  it  is  shown  that  the  method 
converges  about  ten  times  faster  than  the  steepest  descent 
procedure  mentioned  above.   The  new  calculations  of  the 
contraction  coefficient  produced  by  the  reflection  program 
are  also  presented  in  this  section.  They  agree  extremely  well 
with  those  of  Bloch  [1]  . 
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The  evidence  offered  herein  clearly  demonstrates 
the  superiority  of  this  scheme  in  handling  the  vena  contracta 
model,  and  we  are  led  to  conclude  that  reflection  is  the 
appropriate  method  for  solving  a  wider  class  of  free 
boundary  problems  which  arise  in  fluid  dynamics  and  plasma 
physics.   Furthermore,  the  results  tend  to  encourage 
the  use  of  reflection  rules  whenever  they  are  applicable 
in  solving  more  general  nonlinear  elliptic  boundary  value 
problems  in  two  independent  variables . 

As  a  final  note,  we  present  in  Appendix  I  a 
re-evaluation  of  a  calculation  of  the  axisymmetric  contrac- 
tion coefficient  by  a  perturbational  technique  [3] ,  which 
motivated  the  present  study. 
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Part  2.   REFLECTION 

A.   Formulation  of  Free  Surface  Problems 

One  of  tne  principal  sources  of  nonlinear  boundary 
value  problems  for  which  the  reflection  scheme  is  suitable 
is  the  free  surface  problem  in  fluid  dynamics ,  so  a  survey 
of  the  known  results  in  this  area  will  serve  as  our 
introduction  to  the   method. 

One  must  solve  a  partial  differential  equation  for 
the  stream  function  ip  in  a  region,  called  the  flow  region, 
occupied  by  the  fluid.   This  region  is  bounded  partly  by 
fixed  walls  and  obstacles  and  partly  by  a  constant-pressure 
medium,  such  as  the  atmosphere.   The  locations  of  the  walls 
and  obstacles  are  known  data,  but  the  location  of  the  free 
surface  must  be  determined  as  part  of  the  solution.   The 
value  of  ip  is  given  on  all  boundaries,  and  some  other 
condition  (usually  expressing  the  fact  that  the  pressure 
in  the  fluid  must  match  the  atmospheric  pressure)  is  given 
on  tne  free  surface.   Thus  one  is  presented  with  a  differential 
equation  to  solve  in  an  unknown  region  with,  usually,  nonlinear 
boundary  conditions. 

Specifically,  in  two-dimensional  potential  flow 
problems,  the  components  of  velocity  in  the  (x,y)  plane 
are  obtained  from  i>    through  the  equations 

V   =  rr^        ,  V   =  -  -K^ 

X    dy  y       dX 
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The  equation   for  ij;  then  expresses  the  condition  that  the 
flow  is  irrotational: 

-  curl  V  =  A4;  =  0  . 

The  flow  on  the  fixed  and  free  boundaries  is  parallel 
to  the  boundaries,  so 

Tp   =    constant 

there.   On  a  free  boundary,  the  pressure  p  is  constant. 
Since  Bernoulli's  law  relates  pressure  and  velocity  [2], 

p  +  p-  V  =  constant  , 

the  velocity  is  also  constant  there.   Thus  on  the  free 
boundaries  where  the  velocity  is  given  by  the  normal 
derivative  3i|i/3n,   we  must  have 

r.  '   =   constant  . 
3n 

Of  course  the  values  of  all  these  constants  must  be  chosen 
to  satisfy  conditions  at  the  boundaries,  at  infinity,  etc. 
In  three-dimensional  problems  with  axial  symmetry, 
we  let  x  denote  distance  along  the  axis  of  symmetry  and 
y  denote  radial  distance  from  this  axis.   The  velocity  is 
derived  from  the  stream  function  i)   according  to  the  equations 

V  =  i-|l   ,        V   =  -  iM  . 
X    y  dy  y      y  dx 

If  the  flow  is  irrotational. 
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curl  v=--(Ai|j--|^)  =0  , 

y      y  °y 


and  therefore 

A4;  -  ^  |i  =  0  . 
y  9y 

Again,  the  flow  on  the  free  and  fixed  boundaries  is 
tangential,  so  on  these  Doundaries 

ip   =  constant  . 

Bernoulli's  law  dictates  constancy  of  speed  on  the  free 
boundaries  as  before.   This  is  now  expressed  by 


1  d\b 

—   T-"-  =  constant  . 

y  dn 


It  has  been  shown  [5,7]  that  these  problems  have 
unique  solutions ,  and  furthermore  that  the  free  boundaries 
are  analytic.   This  latter  fact  will  be  helpful  in 
implementing  the  reflection  scheme. 

An  example  of  this  kind  of  problem,  which  will  be 
treated  in  our  paper,  is  the  vena  contracta.   Water  is 
kept  under  pressure  in  an  infinite  tank  bounded  by  the 
y-axis ,  and  allowed  to  escape  through  an  opening  in  the 
wall.   The  opening  is  a  slit  in  the  two-dimensional  case, 
and  a  hole  in  the  three-dimensional  axisymmetric  case. 
If  we  neglect  gravity,  the  water  escapes  through  the 
aperture  in  a  stream,  which  eventually  looks  like  a  tube 
of  fluid  in  uniform  motion  (see  Fig.  1) .   The  boundary 


-8- 


of  the  stream  is,  of  course,  unknown,  but  the  pressure 
at  this  boundary  must  equal  atmospheric  pressure.  This 
condition  determines  a  unique  solution.   As  we  mentioned 
above,  the  boundary  is  an  analytic  curve  in  the  (x,y)  plane, 

The  two-dimensional  problem  can  be  solved  by  the 
hodograph  method.   The  fiinction  ^p   is  harmonic  and  can  be 
combined  with  its  harmonic  conjugate  c})  to  give  an  analytic 
fxinction  r,   =   (i>   +   i^    .      The  velocity  is  then  given  by 


_   d? 

V 

-i   V 

/ 

X 

y 

dz 

hence   v  -iv   is  analytic  also.   The  regions  in  which 
X   y 

the  velocity  and  c.  are  defined  are  easily  determined, 
and  the  Schwartz-Caristof f el  formula  yields  the  mapping 
between  them,  so  an  integration  of  the  above  equation 
gives  the  function   C(z).   The  result  is 

-(x+iy)  =  -  2i  +  c  +  I  e   2    -  f  (1  +  e   ^   ) ^^^ 

+  2  log  [1  +  (l+e-^^)^/2j 

IT 

for  some  choice  of  the  constants. 

An  interesting  parameter  in  this  problem  is  the 
contraction  coefficient,  which  is  the  ratio  of  the  cross- 
section  area  of  the  stream  at  infinity  to  the  area  of 

IT 

the  aperture.   The  above  equation  yields  C^  =   :j^-j:2"  ^  .611015 
for  the  two-dimensional  case.   The  axially  symmetric  value 
must  be  evaluated  by  approximate  techniques;  in  [1],  Bloch 
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calculated  it  to  be  .59135  +  .00004.   We  will  present 
another  evaluation  below. 

B.   Parametrization  by  a  Conformal  Transformation 

The  solution  of  free  boundary  problems  can  often 
be  aided  by  a  parametrization  via  a  conformal  map.   One 
is  trying  to  solve  an  equation  for  \p    in  terms  of  x  and  y 
in  an  unknown  region.   Now  if  we  regard  the  region  as  a 
domain  in  the  complex  plane  of  z  =  x  +  iy,  and  if  it  is 
simply  connected  and  not  the  whole  plane,  then  the  Riemann 
mapping  theorem  tells  us  that  this  domain  can  be  regarded 
as  the  image,  under  an  analytic  function  z  =  z (w) ,  of 
something  specific,  like  a  rectangle,  in  the  w-plane. 
Furthermore  we  can  specify  the  boundary  correspondences 
of  three  points.   If  we  have  this  map,  we  can  parametrize 
the  equation  for  ij;(x,y)  to  yield  an  equation  for  i|;  as  a 
function  of  u  and  v,  where  w  =  u+  iv;  this  equation  now 
has  to  be  solved  in  a  simple,  known  domain. 

The  parametrization  of  the  ijj  equation  in  the  case 
discussed  above  is  aided  by  certain  identities  involving 
the  gradient  operator  V  and  the  Laplacian  A;  these 
identities  are  immediate  consequences  of  the  Cauchy-Riemann 
equations  for  the  analytic  fxinction  z  (w)  .   If  V   denotes 
the  gradient  operator  with  respect  to  x  and  y  while  V 
denotes  this  operator  with  respect  to  u  and  v,  then  for 
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any  real  twice-dif ferentiable  functions  f(x,y)  and  g(x,y), 
we  have 

V^  f  (x(u,v)  ,y  (u,v)  )  •  V^g(x(u,v)  ,y  (u,v)  ) 

and 

A^  f(x(u,v),y(u,v))   =   1^1   A^  f(x,y)  . 

Therefore,  for  two-dimensional  potential  flow 
problems,  the  equation 


becomes 


A^  '\>U,y)    =   0 


A   <J;(u,v)  =  0  . 
w  ^   ' 


The  boundary  values  of  ^    just  become  transferred  to  the 

images  of  those  boundaries.   To  facilitate  studying  the 

constant-speed  condition,  let  us  introduce  the  symbols 

n  ,  t  ,  n  ,  and  t   to  denote  normal  and  tangential 
z'   z   w       w  ^ 

directions  in  the  z  and  w  planes.   Then  by  the  chain  rule, 

9i|j  _  9  ij;   w    9  j)        w 

8n     9n    9n     9t    9n 
z     w    z     w     z 

But      9ip/9t      =    0   on   the   pre-image   of    the    free   boundary 
because   i|j    is    constant   there.      So    the    constant  speed 
condition  becomes 

9iij  9  \1;  w  .        . 

V-'—-  =  -"r—-  TT =    constant. 

9n  dn  9n 
z             w        z 
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since  z (w)  is  conformal  and  maps  boundaries  to  boundaries, 
the  normal  to  the  bounaary  is  preserved  and  we  can  write 

3n        9n  ^ 


3n      '^  3n        ^  '  3n  ' 
z         w  w 


provided  the  positive  directions  along  the  normals  are 
chosen  appropriately;  thus  we  can  rewrite  our  condition  as 


^r^  /    l-:r^— I  =  constant 

Art         '        i  Ar\      ' 


3n     '  9n 
w      w 

on  the  free  boundary. 

The  axially  symmetric  equation  for  i|j  is  parametrized 

just  as  easily  when  it  is  written  in  the  form 

V  y  '  V  4^ 

AW; 5 ±_  =  0  . 

z^       y 


It  transforms  to 


V  y  •  V  i^ 

A     ]h     -     =  0 

w^       y 


in  the  w-plane .   The  boundary  values  of  ip  are  carried  over 
as  before,  ana  the  constant  speed  condition  becomes 

—  — -  /   TT =  constant. 

y  3n     ' 3n  ' 
-^    w      w 

So  now  the  problem  for  ip   is  simpler,  given  the 
analytic  map;  one  has  an  equation  with  boundary  conditions 
similar  in  form  to  the  original,  but  these  are  to  be 
solved  in  a  known  region. 

The  fly  in  the  ointment  here  is  that  one  must  also 
find  the  conformal  map.   And  when  there  is  an  unknown  boundary. 
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as  in  free  surface  problems,  the  mapping  problem  is 
coupled  with  the  problem  for  ip ,    through  the  free  boundary 
condition  (as  in  \^\    /    \^\    =   1);  so  one  cannot,  in 
general,  solve  the  equation  and  find  the  mapping  separately. 
Therefore  we  shall  examine  the  conformal  mapping  problem 
in  detail. 

C.   Calculation  of  Conformal  Maps 

Let  z  (w)  be  analytic;  then,  of  course,  x  and  y 
are  harmonic  functions  of  u  and  v: 

Ax  =  0  ,      Ay  =  0  . 

Now  reversing  our  point  of  view,  if  we  wish  to  find  a 
conformal  map,  we  can  start  by  regarding  these  as 
differential  equations  to  be  solved  for  x  and  y; 
but  then  the  boundary  conditions  must  be  determined  to 
fix  the  solution. 

First  of  all,  the  image  z  =  x  +i  y  of  a  point  on 
the  boundary  must  lie  on  the  curve  bounding  the  image 
region.   This  is  a  rather  subtle  requirement,  however, 
because  it  does  not  specify  x  and  y  individually;  at  most 
it   gives  one  in  terms  of  the  other,  and  in  the  case  of 
an  unknown  free  boundary,  it  tells  us  nothing.  More 
conditions  are  needed. 

In  [1],  Bloch  utilized  the  theory  of  Plateau's 
problem  to  get  the  necessary  boundary  conditions.   As 
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described  in  [2] ,  if  one  has  continuously  dif ferentiable 
functions  x(u,v)  and  y(u,v)  defined  on  a  rectangle  in  the 
(u,v)  plane  and  mapping  this  rectangle  onto  a  domain  bounded 
by  a  closed  curve,  and  if  one  minimizes  the  Dirichlet- 
Douglas  integral 

I   ((Vx)^  +  (Vy)^)  du  dv 
R 

subject  to  a  three-point  normalization  (i.e.,  the  images 

of  three  points  on  the  boundary  of  the  rectangle  are  fixed) , 

then  the  resulting  function  x+iy  of  the  complex  variable 

u+iv  is  unique  and  analytic.   To  minimize  the  integral 

under  any  auxiliary  conditions,  Bloch  shows,  one  must 

have  Ax  =  Ay  =  0 ,  and  also  on  the  boundary  of  the  rectangle 

the  normal  derivative  of  the  vector  (x,y)  must  be 

perpendicular  to  its  tangential  derivative.   One  can  express 

this  condition  as 

x^x   +  y.y   =  0 
t  n    -^  t-'  n 

where  n  and  t  denote  normal  and  tangential  derivatives 

respectively.   The  left-hand  side  is  the  dot  product  of 

the  vectors  (x  ,y  )  and  (x  ,y^) .   In  the  case  where  the 

rectangle  in  the  u,v-plane  has  sides  parallel  to  the  axes, 

the  orthogonality  condition  asserts  that 

XX  +  y  y  =0 
u  V    -^  u-'  V 

on  all  sides. 
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To  understand  exactly  what  this    last  condition 
says  about  the  solutions  to  the  equations  Ax  =  Ay  =  0 , 
one  must  investigate  it  further.   Following  Bloch,  we 
note  that  the  function  x  x  +y  y   is  harmonic  if  x  and  y 
are,  so  if  it  vanishes  on  the  boundary  then  it  vanishes 
inside  the  rectangle  as  well.   Its  harmonic  conjugate  is 

1,2^2     2     2.    Ij, 
^(x   +y   -X  -   Y    )    =    -^  I 
2    V    •'v     u    ■'u     2 

as  direct  calculation  shows ,  so  this  latter  function  must 
be  constant  in  the  rectangle.  Furthermore,  if  this 
constant  is  zero,  then  the  map  x+  iy  is  conformal 
because  it  satisfies  the  Cauchy-Riemann  equations. 

In  summary,  if  the  Dirichlet-Douglas  integral  is 
minimized  subject  to  a  three-point  constraint,  the  map  is 
analytic  and  the  constant  His  zero;  and  if  it  is  minimized 
subject  to  additional  constraints  but  with  H  =  0 ,  the 
map  is  still  analytic. 

Bloch  actually  solves  the  problem  by  enforcing 
a  four-point  condition  specifying  the  images  of  the  corners 
of  the  rectangle  and  iterating  his  answer  by  a  steepest 
descent  correction  according  to  the  following: 

1.  moving  (x,y)  along  the  boundary  in  such  a  way 

as  to  drive  x  x  +y  y   to  zero, 
u  V  -^  u-*  V 

2.  moving  the  free  boundary  itself  in  such  a  way 
as  to  drive  the  pressure  to  its  correct  value, 

3.  moving  the  (prescribed)  position  of  two  of  the 
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corners  in  such  a  way  as  to  drive  H  to  zero. 
This  scheme  is  successful,  but  it  is  our  ooal  here  to 
establish  that  reflection  is  equally  successful  and  much 
more  efficient. 

D.   Reflection 

To  describe  the  use  of  reflection  in  seeking  a 
conformal  map  numerically,  let  us  consider  a  rectangle 
in  the  w-plane.   It  is  to  be  mapped  onto  a  region  in  the 
z-plane  bounded  by  real  analytic  curves  (Fig.  2) .   To 
find  the  map,  we  observe  that  x  and  y  are  harmonic  functions 
of  u  and  v,  and  we  use  this  in  a  finite  difference  scheme. 

Let  us  first  set  up  a  mesh  in  the  rectangle  (Fig.  3) . 
If  we  identify  the  points  by  integer  pairs  (i,j),  enumerating 
them  in  the  u  and  v  directions   respectively,  then  for 
harmonic  functions  we  have  the  approximation 

1  2 

X.  .  =  4  (x.,,  .+x.  ,  .+x.  .,,+x.  .  ,)  +  0(h  ) 
1,3    T   1+1, D   1-1, D   i,D+l   i,D-l 

where  h  is  the  mesh  size,  assumed  uniform.   Of  course, 
a  similar  equation  holds  for  y.   This  is  the  difference 
equation  approximating  the  equation  Ax  =  0 . 

This  finite  difference  analogue  of  the  underlying 
differential  equation  can  be  written  at  every  interior  mesh 
point  of  our  rectangular  grid.   In  standard  treatments  of 
boundary  value  problems,  e.g.  the  Dirichlet  problem,  another 
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difference  equation  approximating  the  boundary  condition 
is  written  at  the  mesh  points  on  the  edge  of  the  rectangle. 
However,  in  nonlinear  examples  such  as  the  conformal  mapping 
of  the  rectangle  onto  some  arbitrary  region  this  may  not 
be  feasible.   To  circumvent  the  difficulties  we  introduce 
new  mesh  points  just  outside  the  rectangle  and  use  them 
to  write  this  same  difference  analogue  of  the  differential 
equation  at  the  boundary,  as  well  as  the  interior,  mesh 
points.   Then  the  problem  becomes  finding  reflection 
rules,  based  on  the  boundary  conditions,  which  will 
determine  the  values  of  the  solution  at  the  reflected 
mesh  points.   This  procedure  is  familiar  in  treating  the 
Neumann  problem,  and  it  usually  leads  to  a  convergent 
iterative  scheme  for  finding  the  solution  even  when 
nonlinearities  occur. 

Clearly  the  technique  is  applicable,  in  theory 
at  least,  to  any  problem  in  which  the  solution  can  be 
continued  across  the  boundary  according  to  known  rules. 
In  [8] ,  Lewy  demonstrates  this  property  for  a  class  of 
analytic  elliptic  equations,  using  the  appropriate  Riemann 
functions.   In  our  case,  where  a  straight  line  in  the 
w-plane  is  mapped  onto  an  analytic  curve  in  the  z -plane, 
the  property  of  analytic  continuation  is  a  well  known 
result  of  the  theory  of  conformal  mapping  [9] . 

So  now  we  form  additional  mesh  points  outside  the 
rectangle  which  are  reflections,  through  the  sides,  of 
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those  mesh  points  adjacent  to  the  sides,  and  x+iy  can  be 
defined  at  these  reflected  points  (at  least,  for  small 
enough  h)  (Fig.  2),   Thus  we  can  write  the  numerical  analogue 
for  Laplace's  equation  at  each  boundary  point,  even  at  the 
corners;  each  one  has  a  left  and  right  and  an  upper  and 
lower  neighbor. 

This  leaves  us  with  the  problem  of  writing 
equations  for  the  values  of  x  and  y  at  the  reflected  points, 
that  is,  the  formula  for  the  reflection  rule.   To  derive  this 
rule,  we  must  consider  more  carefully  the  exact  description 
of  the  analytic  curves  bounding  the  image  region  in  the 
z-plane . 

Let  us  say  the  real-analytic  curve  T    is  one  of 
the  bounding  curves  in  the  z-plane,  and  its  pre-image  is 
the  bottom  of  the  rectangle,  v  =  0  (Fig.  2).   The  equation 
for  r  is  of  the  form 

g(x,y)  =  0  , 

where  g  can  be  expanded  locally  in  a  real  power  series. 
Making  the  substitutions 

z+z  z-z 

(where  a  bar  denotes  complex  conjugation)  we  get  the 
equation  for  T    in  the  form 

F(z,z)  =  0  , 

where  F(z, ,z„)  is  complex-analytic  in  z,  and  z_  because 

F  has  a  power  series  expansion,  namely,  that  derived  from  g. 
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Now  if  w_  is  an  interior  point  in  the  w  plane  and  w   is  the 
reflected  image  of  w_,  reflected  through  the  line  v  =  0,  and 
if  z_  =  z(w_)  and  z   =  z (w  ) ,  then  we  assert  that 

F(z^,z_)  =  0  . 

This  is  the  basic  formula  which  gives  values  of  z  at  the 
reflected  points  in  terms  of  the  values  at  interior  points. 
Before  we  prove  this,  let  us  examine  two  special 
cases . 

Case  1.   r  is  a  straight  line,  say  x  =  1. 
Rewriting  this  as  x-1  =  0,  then  as 

we  get 

z  +z 
F(z^,Z2)  =  -V^-  1 

We  know  the  rule  for  reflecting  in  this  case;  it 


is  (z  -1)  =  -(z_-l).   This  agrees  with  F(z, ,z„)  =  0 

2     2 
Case  2.   F  is  a  circle,  say  x   +  y   =1. 

This  is  written  as  zz-1  =  0  and 

F(Zj^,Z2)  =  z^Z2  -  1  . 

The  rule  for  reflection  is 

1 
z_ 

again  agreeing  with  F(z  ,z_)  =  0. 
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Of  course,  if  F(z,z)  =  0  is  to  specify  a  curve, 

then  Re  F(z,z)  and  Im  F(z,z)  are  not  independent  functions; 

if  they  were,  the  (complex)  condition  F(z,z)  =  0  would 

define,  as  a  locus,  the  intersection  of  two  different 

curve:;,  i.e.,  certain  isolated  points.   In  the  two  cases 

considered  above,  Im  F(z,z)  was  identically  zero,  as  it 

will  always  be  if  F  is  derived  from  g  in  the  manner 

described. 

Now  we  prove  the  rule  for  reflection.   Define 
z  and  w  as  before. 

Theorem.   If  z (w)  is  analytic  and  if  it  is  defined  on  both 
sides  of  the  line  v  =  0,  mapping  this  line  onto  the  analytic 
curve  F(z,z)  =  0,  then  F(z  ,z_)  -  0. 

Proof:   Noticing  that  w   and  w_  are  related  by  w_  =  w  , 
we  define  the  function 


f(w)  =  F(z(w)  ,zCw)) 


where  z (w)  denotes  the  conjugate  of  z (w) .  When  w  =  w  , 
w  =  w_  and  f (w)  is  just  F(z  , z_) .  We  know  F  is  analytic 
in  its  first  and  second  variables.   Since  z (w)  is  also 


analytic,  z (w)  is  analytic  and  thus  f  is  analytic  in  w, 

by  composition.   When  w  lies  on  the  bottom  of  the  rectangle, 

w  =  u  =  w,  then  z  lies  on  T,    so 

f(w)  =  F(z{u),i(u))  =  F(z,z)  =  0  . 

Therefore  f (w) ,  an  analytic  function  of  w,  vanishes  on 
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the  line  v  =  0;  hence  it  is  identically  zero.   Thus 
F(z^,z_)  =  f(w^)  =  0  .  Q.E.D. 

From  the  above  proof,  we  can  readily  see  that  the 
same  rule  holds  for  reflection  through  the  other  sides  of  the 
rectangle  also,  i.e.  F(z^,z_)  =  0  is  the  general  rule  for 
analytic  continuation . 

Now  we  have  a  feasible  scheme  for  calculating 
conformal  maps.   We  write  the  numerical  equivalent  of 
Laplace's  equation  at  each  interior  and  boundary  point, 
and  we  write  F(z  , z_)  =  0  for  each  reflected  point. 
The  next  section  discusses  the  consequence 
procedure . 

E.   Implications  of  the  Reflection  Rules 

Now  we  must  reverse  our  point  of  view  in  this  sense; 
we  set  up  Baplace's  difference  equation  at  each  interior 
and  boiindary  point  of  the  rectangle,  and  we  write  the 
reflection  rule  for  each  of  the  reflected  points .   What 
can  we     say  about  the  solution  to  this  system? 

For  the  moment  let  us  assume  that  the  solution 
exists  and  is  unique  for  every  sufficiently  small  mesh  size  h, 
and  that  the  solution  converges  to  a  function  z(u,v)  =  x+iy, 
where  x  and  y  are  harmonic,  as  h  goes  to  zero.   These 
assumptions  will  be  discussed  in  the  next  section.   What  are 
the  properties  of  s  on  the  boundary? 
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Let  us  examine  a  particular  point  w„  on  the  boundary 
of  the  rectangle,  together  with  its  inner  neighbor  w_  and 
its  reflected  neighbor  w   (see  Fig.  2) .   As  h  goes  to  zero, 
the  three  points  z   =  z  (w  )  ,  z-  =  z(Wq),  z_  =  z  (w_)  ,  must 
all  converge  to  the  same  point  z, .   But  F(z  ,z_)  =  0,  so 
F(z,,z^)  =  0  and  z,  lies  on  F.   Another  way  cf  seeing  this 
is  to  notice  that  z   and  z_  always  lie  on  opposite  sides 
of  r,  so  z   must  lie  on  F .   We  conclude  that  the  limiting 
function  z  maps  the  boundary  of  the  rectangle  onto  the 
boundary  of  the  region  in  the  z-plane. 

Now  it  is  our  goal  to  prove  that  the  normal  and 
tangential  derivatives  of  z  are  perpendicular,  i.e., 
X  X  +y  y  =  0.   To  this  end  we  must  study  the  geometry  of 
the  curve  T . 

First  v^e  assume  that  the  function  F(z,,Z2)  was 
derived  from  a  real  function  g(x,y)  as  described  earlier, 

so  that  F(z,z)  is  always  real.   If  z  is  given  as  a  function 

dF 
of  a  real  variable  s,  z  =  z(s)  ,  then  -a—  must  be  real. 

Therefore 

dF 

5^  =  F,z   +  F-z 

ds     Is     2  s 

is  real.   Since  z(s)  is  arbitrary,   we  conclude 
Re  F,  =  Re  F-,  Im  F,  =  -  Im  F_;  in  other  words 


Now  if  z(t)  is  a  parametrization  of  the  curve  T, 
F(z(t)  ,z(t)  )  =  0.   Thus 


-22- 


g  =  0  =  F^z^+  F^i^  =  F^z^^  F^i^  =  2  ReCF^z^) 


0  =  Re  F^x  +  Im  F2y^ 

Now  (x  ,y.)  is  tangent  to  the  curve  r,  so  (Re  Y^,    Im  F2)  is 
normal  to  the  curve  and  (-Im  F2,Re  F2)  is  again  tangent  to  r. 

If  we  associate,  with  each  complex  number,  a  vector 
v/hose  X  and  y  components  are  the  real  and  imaginary  parts 
respectively,  we  can  say  that  F2(^f^)  is  normal  to  the 
curve  F(?,t)  =  0. 

We  can  get  more  information  about  r  by  applying 
the  analytic  form  of  the  implicit  function  theorem.   If 
F(z,,Z2)  is  analytic  in  both  variables  with  F(a,b)  =  0 
and  F,(a,b)  ^   0,  then  F(z,,Z2)  =  0  defines  z^  as  an  analytic 
function  of  z„  in  a  neighborhood  of  (a,b) ; 

2 
(z,-a)  =  Cq  +  c,(z2-b)  +  C2(z2-b)   +  ...  . 

^2 
Furthermore,  c„  =  0  and  c,  =  -  ^  (a,b) . 

We  have  F(z  , z  )  =0  defining  z   in  terms  of  z_. 

Let  us  suppose  that  z   and  z_  both  approach  ?  as  h  -»■  0; 

we  just  proved  r>   was  on  the  curve  T;  so  F(^,C)  =  0. 

The  theoren  then  tells  us  that  we  can  write 

F2(?,C)  -   -  2 

(z.-s)  =  -  — (z  -I)    +  c^(z_-C)   + 

F  (CO 


•  •  • 


Using  F,(5,5)  =  Y^it.iO    and  adding  ^-z_  to  both  sides, 
we  ultimately  have 
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F  _ 


As  the  mesh  size  h  of  the  rectangle  becomes  smaller, 

?-z_       ^+~^_ 
the  quantities  -,- and   — -r- should  give  approximations 


to  the  derivative  of  z  in  the  direction  norraal  to  the  side 

of  the  rectangle,  denoted  z  .   Thus  if  we  divide  the  above 
-^  n 

equations  by  2h  and  take  the  limit,  we  get 

^2       -   - 

z   =  T  [F^z  +F~z  ] 

"    i]^    2  n   2  n 

Since  the  coefficient  of  F^  is  real,  the  direction  of  the 

vector  z  ,  if  it  is  non-zero,  is  the  same  as  that  of  F-, 
n'  '  2 ' 

i.e.   normal  to  the  curve  T.       (A  case  where  z   is  zero 

n 

will  be  mentioned  in  the  next  section.) 

Of  course,  the  derivative  z   of  z  along  the  side 
of  the  rectangle  will  be  parallel  to  the  tangent  to  T, 
so  we  will  have 

^n^t  -^  ^n^t  =  ° 

for  the  function  z(u,v)  along  the  boundary  of  the  rectangle. 
Reasoning  as  before  for  these  harmonic  functions  x  and  y, 

we  know  that  the  latter  equation  holds  throughout  the 

2     2     2     2 
rectangle  and  that  x   +y   -x   -y   =His  constant  there. 
^  V   ■'v    u   -^u 

From  what  v/e  have  said  so  far  we  cannot  conclude 
that  z  is  analytic,  i.e.,  that  H  =  0.   For  example,  the 
mapping  x  =  2u,  y  =  v  satisfies  all  the  above  requirements 
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and  is  not  analytic.   To  get  an  analytic  map  we  must  impose 
some  additional  condition  to  drive  H  to  zero;  hopefully 
the  particular  application  will  suggest  such  a  condition. 

F.   Existence,  Uniqueness,  and  Convergence 

We  shall  now  discuss  the  existence  and  uniqueness 
of  the  solution,  and  its  convergence  as  h  -*■  0 ,  on  the 
basis  of  nmmerical  experiments. 

The  method  has  been  programmed  for  several  different 
cases  of  image  regions  bounded  by  four  analytic  curves. 
The  vena  contracta  described  below  leads  to  just  such  a 
region.   The  reflection  rules  were  applied  so  as  to  map 
each  side  of  the  rectangle  onto  one  of  the  bounding  curves. 
This  amounts  to  a  four  point  normalization  because  each 
corner  must  be  carried  into  the  intersections  of  these 
curves.   The  equations  were  iterated  using  overrelaxation 
to  accelerate  convergence  and  in  each  case  the  iterates  did 
converge.   This  indicates  (but  does  not  prove)  that  the 
solution  of  the  four  point  normalization  problem  exists 
and  is  unique,  for  if  there  were  two  solutions,  one  would 
expect  the  iterates  to  oscillate  between  them.   By  shifting 
one  of  the  curves  one  can  get  different  values  for  H; 
only  one  position  makes  H  =  0,  and  then  the  finite  difference 
analogues  of  the  Cauchy-Riemann  equations  are  satisfied. 

No  proof  of  convergence  of  the  method  as  h  ^  0 
has  been  developed,  but  the  problem  has  similarities  to  the 
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Dirichlet  and  Neijtmann  problems,  whose  finite  difference 
analogues  have  been  analyzed  [2,6,11].   In  the  latter  cases, 

the  solutions  of  the  difference  equations  approach  those 

2 
of  the  differential  equaton  with,  an  error  of  the  order  0(h), 

so  we  expect  the  same  of  our  scheme  also.   Nimerical 

experiments  confirm  that  the  convergence  is  0(h  )  (see  below) . 

Let  us  briefly  examine  the  three  point  normalization. 

If  we  choose  to  map  a  rectangle  onto  a  domain  bounded  by 

three  analytic  curves,  we  would  naturally  map  two  adjacent 

sides  onto  one  of  the  curves,  and  we  would  not  be  specifying 

the  location  of  the  included  corner.   In  Figure  4  we  have 

labeled  the  corner  w  ,  while  its  neighboring  mesh  points  on 

the  top  edge  are  denoted  w,  and  w^,  and  on  the  vertical  edge, 

w^  and  w..   The  neighboring  interior  point  is  called  w_,  and 

its  reflected  images  through  the  top  and  side  are  called 

w f  '  and  w|   ,  respectively.   Observe  that  the  two  reflections 

of  w_  are  mapped  into  one  and  the  same  point,  z  ,  by  the 

rule  F(z  , z(w_) )  =  0,  because  the  same  reflection  rule  is 

applied  to  both  sides  of  the  rectangle.   If  we  denote  the 

image  of  each  subscripted  w  by  a  z  with  the  same  subscript, 

we  write  Laplace's  difference  equations 

^1=   i   (z_+Zc+2++23) 
so  4z,-z^  =  4z2-z..   Therefore  we  can  write 
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3z  -4z,+z^    3z  -4z„+z. 
c    1   3  _    c    2   4 

2H 2H • 

The  left  and  right  hand  sides  are  finite  difference 

2      3  z      8  z 
approximations  of  order  0(h)  to  ^—  and  ^ —  at  the  corner, 

so  if  we  can  assume  that  these  converge  to  the  derivatives, 

we  have  z  =  z   there.   Thus  x  =  x  and  y  =  y  ,  so 

u    V  u    V     -'u   -^v' 

2   2   2   2 
H  =  X  +y  -X  -y   is  zero  at  the  corner,  hence  zero  evervwhere, 
V  ■'v  u  ■'u 

and  therefore   z  is  analytic.   We  conclude  that  if  the 
reflection  scheme  converges  for  the  three  point  normalized 
mapping  on  a  rectangle,  the  limiting  function  is  analytic. 

Notice  that  nowhere  in  the  above  proofs  do  we  assume 
that  the  image  of  boundary  mesh  points  lie  on  the  boundary  of 
the  image  region.   Only  in  the  limit  as  h  -^  0  can  we  be  sure 
of  this  correspondence  of  boiuxdaries . 


G.   Summary 

The  results  of  the  preceding  section  can  be  stated 
concisely  as  follows. 

Consider  analytic  boundary  curves  described  by 
equations  of  the  form 

F(z,z)  =  0  . 

We  write  Laplace's  difference  equation  for  x  and  y  in  a 
region  in  the  w-plane;  the  difference  equations  are  vritten 
for  each  interior  and  boundary  point.   Then  we  write 


F(z^,z_)  =  0 
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for  tlie  reflected  points.   This  procedure  constitutes  the 
reflection  scheme  for  conformal  mapping  problems. 

When  the  w-region  is  a  rectangle  and  the  image  is 
bounded  by  four  analytic  curves,  each  corresponding  to  a 
side  of  the  rectangle,  our  numerical  experiments  lead  us 
to  conjecture  the  following  theorem. 

Theorem.   The  solutions  x,  and  y,  of  the  equations  of  the 
reflection  scheme  converge  to  harmonic  functions  x(u,v) 

and  y(u,v)  as  the  mesh  size  h  is  diminished.   The  rate 

2 
of  convergence  is  0(h  ) . 

Furthermore,  whenever  the  solutions  and  their 

difference  quotients  converge,  we  have  proved  the  following. 

Theorem.   Under  the  convergence  assumptions  stated,  the 
function  z (w)  maps  the  sides  of  the  rectangle  onto  the 
boundary  curves,  and  on  these  sides  the  derivatives  of  z 
in  the  normal  and  tangential  directions  are  perpendicular 
to  each  other. 

As  we  mentioned,  it  follows  immediately  that 

XX  +  y  y  =0 

u  v   -^  u-'  V 

in  the  rectangle,  and  therefore 

2     2     2     2 
X  +  y   -  x   -  y   =  constant  . 
V    -'  V     u     u 

Experimentation  shows  that  the  value  of  this 

constant  depends  on  the  locations  of  the  curves  in  the 

z^plane,  and  the  constant  can  be  made  zero  by  choosing 
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the  locations  appropriately.   This  condition  will  then 
make  z (w)  analytic. 

The  above  procedure  seems  to  be  most  suitable  for 
applications  of  the  method  of  reflection  when  the  curves 
in  the  z-plane  are  analytic  and  known.   However,  as  we 
proceed  in  Part  3  it  will  become  clear  that  neither  of 
these  conditions  is  essential,  and  one  can  often  adopt 
the  technique  to  less  restrictive  cases. 
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Part  3.   THE  VENA  CONTRACTA 

The  reflection  technique  will  now  be  employed  to 
solve  the  two-  and  three-dimensional  vena  contracta  problems. 
This  section  describes  the  details  of  the  procedure. 

A.   The  Physical  Problem  and  Parameters 

As  described  earlier,  for  the  physical  situation  in 
the  vena  contracta,  one  has  an  infinite  volume  of  incompressible 
fluid  under  pressure  behind  an  infinite  plane  wall;  the  wall 
has  an  aperture  which  is  either  a  slit  (two  dimensions)  or 
a  circular  hole  (three  dimensions) .   A  jet  of  fluid  escapes 
from  the  aperture  and  its  cross-section  contracts  as  it  moves 
away  from  the  wall;  asymptotically,  it  becomes  a  tube  of  fluid 
in  uniform  motion. 

In  eitlier  case,  one  needs  only  to  examine  the  flow 
in  an  (x,y) -plane.   For  the  infinite  slit,  one  has  transla- 
tional  symmetry  in  the  z-direction;  for  the  circular  hole, 
one  has  rotational  symmetry  about  the  x-axis . 

The  motion  is  described  by  the  stream  function  'Jj(x,y), 
which  gives  the  x  and  y  components  of  velocity  by 

V  =   il  /  y^ 

X     3y  '^  -^ 

V  =  -  Mi  /  y"^ 
y     3x  ^  -^ 

where  m  =  0  in  two  dimensions  and  m  =  1  in  three  dimensions. 
The  equation  determining  41  is 

A^  -  E  |i  =   0  . 
y  9y 
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Let  tiie  y-axis  be  the  wall,  with  the  jet  issuing 
to  the  right  and  the  aperture  extending  from  y  =  Yq  to 
y  =  -Y  .   Then  the  motion  is  symmetric  about  the  x-axis 
and  we  can  merely  consider  the  problem  in  the  upper  half 
plane . 

Since  the  x-axis  is  a  streamline,  we  can  set  i|^  =  0 
there  arbitrarily.   The  wall  and  the  free  surface  also  form 
one  streamline,  so  ^p   will  be  a  constant,  iJj-,  there-  The 
fluid  behind  the  v/all  is  under  a  pressure  p^  at  infinity, 
and  atmospheric  pressure  in  front  of  the  wall  is  a  constant, 
p- .   The  free  streamline  approaches  a  horizontal  line,  y  =  Y^, 
asymptotically.   The  speed  of  the  flow  on  the  surface  must 
be  a  constant,  U^ ,  by  Bernoulli's  law,  so  if  n  denotes  the 
direction  normal  to  the  streamline,  we  have 

ii=  u 

9n     0 


in  two  dimensions,  and 


1  3'lJ  _  U 


in  three  dimensions  on  the  free  surface. 

These  constants  are  related,  of  course.   The  speed 
U»  depends,  via  Bernoulli's  law,  on  the  difference  p^-  p 
The  value  i|^^  determines  the  rate  at  which  fluid  escapes  in 
the  jet,  so  it  depends  on  the  other  parameters  also. 
To  obtain  a  consistent  set  of  values,  we  observe  the  following, 
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First  of  all,  the  behavior  of  ^i   at  infinity  in  the 
jet  is  simple.   In  (m+2)  dimensions. 


for  constants  a   and  3.   This  follows  from  the  differential 
equation  and  the  fact  that  ip  will  ultimately  not  depend  on  x. 
Also,  for  large  x 

7n  -  ^  • 
Thus  the  boundary  conditions  imply 

a  =  Uq 
6  =  Uq/2 

Since  ^   approaches   i|^^  when  y  =  Y^  and  x  increases,  we  have 

Uq^co  ,   m  =  0 

^—  .   ni  -  1 

Therefore  ij^_  can  be  calculated  from  Y^  and  (p^-p  )  . 

Secondly,  it  is  clear  on  physical  grounds  that 
specifying  Y_ ,  p^,  and  p.  completely  determines  the  solution; 
these  are  precisely  the  variables  which  the  experimenter  can 
control.   Furthermore,  only  the  difference  p^-p^^  can  be 
relevant;  an  equal  increase  in  pressure  on  both  sides  of 
the  wall  cannot  change  the  flow  of  an  incompressible  fluid. 
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Third,  the  geometry  of  the  flow  region,  in  particular 
the  location  of  the  free  surface,  depends  only  on  Yq.   Given 
one  flow  with  its  associated  Y^  and  P<„-Pq  (determining  Uq  and 
ip-)  ,  we  can  find  another  flow  for  a  different  value  of  the 
pressure  drop   Poo~Pn   ^^"^  with  tlie  same  Yq  ,  by  calculating 

u'  from  (p1-Pq)  ,  then  i|Jq  from  Uq  and  the  unchanged  Y^,  and 

I 

finally  multiplying  i|j  by  the  constant  ^Q/i>Q-      This  flow 

issues  from  the  same  aperture  Y   and  satisfies  the  pressure 
condition  appropriate  to  p^-p^r  ^"^  because  of  the  uniqueness 

00     [J 

theorem  fV],  it  is  the  only  one.   It  utilizes  the  same  flow 

region,  hence  the  same  free  surface.   So,  as  we  claimed, 

Yf.  alone  determines  the  free  surface,  and,  consequently,  Y^. 

Fourth,  notice  that  a  simple  magnification  of  the 
flow  region  for  Y^  gives  a  flow  region  for  a  different  Y- 
(and  hence  the  unique  one)  ,  so  the  non-dimensional  ratio  Y^Yq 
is  independent  of  any  of  the  parameters  of  the  problem.   The 
number  iX^/Y^)^'^^    (m,  as  before,  is  0  or  1  in  two  and  three 
dimensions  respectively)   is  called  the  contraction  coefficient, 
and  we  shall  determine  this  number  in  both  the  two-  and  three- 
dimensional  cases. 

It  follows  that  we  can  determine  the  flow  region  by 

specifying  Y^  instead  of  Y-;  furthermore,  this  convention 

will  facilitate  programming  our  solution,  and  we  adopt  the 

following  normalizations. 

Y  =  1  and  U„  =  1  in  both  two  and  three  dimensions. 

oo  0 

It  follows  that  ^n   =  }    ,-'   "^~,  .   p_  and  p-  are  not  needed 

^0    ]  .  5 ,  m=  1    '^ "°      0 

for  the  solution  (since  we  already  have  U_) ,  and  Y^  must  be 
calculated  as  part  of  the  problem. 
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B .   The  Conf ormal  Map 

As  described  in  Part  1,  the  problen  will  be  attacked 
using  a  confoiT:.al  map.   Theoretically,  we  would  like  to  map  an 
infinite  rectangle  onto  the  flow  region.   The  part  of  the 
rectangle  extending  to  infinity  would  correspond  to  the  jet, 
the  top  of  the  rectangle  would  represent  the  free  boundary, 
and  the  bottom  would  represent  the  x-axis .   The  remaining 
side  would  be  mapped  onto  the  fixed  wall,  and  a  simple  pole 
at  the  lower  corner  would  map  this  corner  to  infinity  in  the 
second  quadrant.   These  correspondences  are  depicted  in 
Figure  5,  where  ABE  is  the  infinite  rectangle. 

In  practice,  we  will  approximate  this  situation 
with  a  finite  rectangle,  truncating  ABE  with  the  vertical 
side  CD.   The  boundary  correspondences  will  now  be  defined 
precisely. 

The  corner  A  of  the  rectangle,  at  the  origin  of  the 
w-plane,  will  be  mapped  to  infinity  in  the  second  quadrant. 
This  requires  a  pole  in  the  transformation  z (w) ,  so  near 
w  =  0  we  have 

z(w)  ^  R/w 

with  a  negative  residue  R  which  must  be  calculated.   The 
side  AB  will  be  mapped  onto  the  wall  from  y  =  <»  to  y  =  Y. . 
The  side  BC  is  mapped  onto  the  free  boundary  r.   So  the 
angle  at  the  corner  B  is  expanded  from  Tr/2  to  it;  this  is  done 
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1/2 
because  near  the  corner  ip  has  an  expansion  in  (z— z„)  '  , 

and  the  doubling  of  the  angle  which  we  will  get  keeps 

2 
the  truncation  error  down  to  0 (h  );  see  [1]  for  details. 

The  side  AD  is  mapped  onto  the  x-axis . 

The  image  of  the  line  CD  is  a  curve  CD'  which 
becomes  asymptotically  a  straight  segment  as  the  rectangle 
is  lengthened.   We  truncate  this  behavior  and  ass\ame  C^D* 
is  a  straight  line,  but  its  location  must  still  be  determined 
(Fig.  5) . 

The  treatment  of  all  the  boundary  conditions  around 
the  rectangle  is  a  complicated  affair,  and  we  will  devote 
the  next  few  sections  to  this  matter. 

C.   Application  of  the  Reflection  Laws  to  the  x-axis 
and  the  Wall. 

The  sides  AD  and  AB  are  each  to  be  mapped  onto  known 
curves,  so  the  reflection  technique  is  appropriate  here. 

AD  is  mapped  to  the  x— axis,  whose  equation  is  z-z  =  0, 
Indicating  images  of  interior  points  adjacent  to  AD  by  x_  and 
y_,  and  their  reflections  by  x  and  y   (Fig.  5) ,  we  derive 

x^  =  x_    and    y^  =  -y_  . 

AB  is  mapped  into  the  y— axis,  and  we  have  a  similar 
situation.   The  reflection  rule  becomes 

y^  =  y_    and    x^  =  -x_  . 
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Arguing  as  in  Part  1,  we  conclude  that  these  equations 
will  demand  that  the  boundaries  do  correspond  and  that 


XX  +  y  y   -  0 
u  V    -^  u-'  V 


along  AD  and  AB , 


D.   The  Reflection  Law  at  the  Free  Surface 

BC  is  to  be  mapped  onto  the  free  boundary  r,  but  the 
situation  is  more  complicated  than  that  treated  in  Part  1. 
The  location  of  the  free  boundary  is  unknown,  and  we  have 
no  simple  equation  like  F(z,z)  =  0  to  work  with.   However, 
we  can  derive  a  r^^flection  rule  from  the  integral  equation 
for  ijj  in  two  dimensions,  and  we  will  modify  it  to  make  it 
suitable  in  three  dimensions.   Then  we  v/ill  have  to  study 
the  implications  of  such  a  procedure  all  over  again. 

To  obtain  the  integral  equation  for  i) ,   we  again 
employ  the  implicit  function  theorem  as  in  Part  2  to  rewrite 
the  equation  F(z,z)  =  0  in  the  form 

z  =  g(z) 

where  g  is  an  analytic  function.   With  this  formulation, 

Garabedian  has  established  in  [3]  that  the  stream  function 

for  gener.ilized  axially  symmetric  flow  in  (m+2)  dimensions 

is  given  by 

z 
1   i  ^:^_^^r"/2.._^,^NNm/2 


4;  (x,y;m)  =  Im  — 
2m 


^1 


(z-t)"'^^(z-g(t)) 

^(_m^_m.  ^.  (z-t)(z-g(t))  )g.(^)l/2  ^^^ 


(z-t)  (z-g(t)) 
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where  z   is  any  point  on  the  free  surface  T,    and  .j~    is 
the  hypergeometric  function 


~  r(i-  ^)^ 

•^   ^   ^        r(-  J)^  j=o  r(j+i)^ 


T         j 


The  integrand  on  the  right-hand  side,  considered  as  a  function 
of  z,  is  essentially  the  Riemann  function  of  the  equation 

8z  3z    2(z-z)        dz 


for 


^(x,y)  =  4^(^7^/  ^^f)  / 


and  this  equation  is  just  the  complex  forr.  of 

Y   W 

which  is  the  usual  equation  for  the  stream  function. 
The  integrand,  considered  as  a  function  of  t,  is  analytic 
so  the  choice  of  tlie  path  of  integration  is  immaterial. 
Finally,  we  add  that  the  normalization  for  ^p    in  this 
equation  differs  from  ours,  where  4/  =  1  on  r,  but  tlie 
additive  constant  will  drop  out  in  our  application. 

The  equation  is  particularly  useful  in  two  dimensions, 
m  =  0 ;  it  becomes 

ij;(x,y)  =  Im  j   g'  (t)     dt  . 

^1 
In  two  dimensions,  i)   is  harmonic.  Thus  it  can  be  combined 
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with  its  hanT\onic  conjugate  (()(x,y)  to  give  an  analytic 

function  Z   =   <t>   +   i   ^ .      The  above  formula  immediately 

displays  this  property,  since  the  right-hand  side  is 

analytic,  and 

z 
<t>   +   i    ^  ^    \      g'(t)^/^  dt  . 

^1 


Differentiating,  we  see  that 

g(z)  = 


z 
C 


(^•)^  dz 


^2 
for  some  z-.   Staying  with  the  two-dimensional  case, 
we  employ  our  usual  notation:   w_,  w„ ,  and  w  denote 
interior,  boundary,  and  reflected  mesh  points  in  the 
rectangle,  while  z_,  z-,  and  z   are  their  images  (Fig.  5). 
The  free  boundary  equation  says 

Zq  -  g(zQ)  =  0 

and  the  reflection  law  F(z  , z_)  =  0  takes  the  form 

z^  -  g(z_)  =0  . 

Siibstituting  for  g(z)  and  subtracting,  we  have  the  reflection 
rule  for  the  two-dimensional  free  boundary  problem. 


^+  -  Zq  =  g(zj  -  g(zQ)  =   j 


z,  -  z-  =  g(z  )  -  g(Zrt)  =   I  ?'   dz  . 


^0 

It  is  now  our  task  to  approximate  this  rule  so  that  it  can 

be  used  in  our  finite  difference  scheme. 
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We  can  integrate  the  right-hand  side  by  employing 

the  power  series  development  of  z,   in  z,  and  that  of  z  in  w. 

I 
Expanding  ?'(z)  about  z  and  using  the  notation  C'(Zq)  =  ^q, 

we  have 


By  the  chain  rule. 


^'  ^  d((j)+ij^l 


dw 


1 
dz/dw 


Zq=z(Wq) 


1   3v      dV    ' 


But  since  ({>   =  ip   =  0  on  the  top  edge  of  the  rectangle. 


^0  =  "^v/'o 


Similarly  we  find 


2 


dz 


'0 


1    r  ^UV      VV 

^      ^0 


'P  z. 
^v  0 

"T77 


where  we  use  the  Cauchy-Riemann  equations  to  eliminate  (j) 
Notice  that  since  ^   is  harmonic. 


^        =   -  ^\, 
^vv     ^uu 


and  the  latter  is,  again,  zero  along  the  top  edge  of  the 
rectangle.  Now  we  use  the  expansion  for  z  in  terms  of  w 
to  derive 
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"2 

dzi        ^-  ^0    ^0^   ^  ^,,2. 

Z^   =   ;;t—       =   -  — , r:—  -  ^  .,     +  0(h   ) 

0    dw'w^       i  h     2ih 

in  terms  of  the  mesh  size  h.   When  all  of  this  is  substituted 

It 

into  the  integral,  the  appearance  of  z.  cancels  and  we  get 

2    2         3 

-h   4i   +ih   iii„tJj  -, 

z^   -   Zq   =   +   0(h    ) 

From   this    approximate   equation  we    formulate   our 

reflect:! on  rule  as 

2  2    3 

-h  ill  +ih  ^   \1) 

\  =  z^  +  — !z !ziHZ 

The  i/^-derivatives  can  be  approximated  by  finite  differences 

accurately  so  that  the  truncation  error  on  the  right  hand 

2 
side  is  of  order  h   higher  th.an  the  left-hand  side  itself, 

which  is,  of  course,  0(h).   Heuristi cally  this  seems  desirable 

because  we  are  using  a  difference  equation  which  approaches 

2 
Laplace's  equation  to  0(h  ) ,  and  we  should  try  to  match  this 

accuracy  in  all  our  other  formulae. 

However,  to  renlly  evaluate  the  effed  of  this  rule, 
we  must  studv  its  consequences  as  we  did  for  the  exact  rules 
in  Part  2.   We  assume  we  have  a  solution  of  Laplace's  differ- 
ence equation  defined  on  the  interior  and  boundary  mesh  points, 
At  the  points  z  which  are  reflected  through  T,  we  have 

2  2    3 

-h  iJj  +ih  li  lii 
-    -       ^v    ^v  uv 
z   —  z_  =  

+     0        ^_~^o 

2 
If  we  multiply  by  (z  -z^^)  ,  divide  by  h  ,  and  take  the  real 

part,  we  get 
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^^  -T E—  =  -  ^  • 

Notice  that  — r— —  nand  — j- —  are  both  difference  approxi- 
mations  to     ^  on  the  boundary,   if  we  assume  convergence 
of  the  scheme  as  we  did  in  Part  2,   we  can  conclude  that 
the  limiting  form  of  the  equation 

_  Tz  9z  ^    ,2 


or  equivalently 


</Nvl   -  1 


is  a  consequence  of  the  reflection  rule. 

There  is  another  condition  imposed  by  this  rule, 
corresponding  to  the  imaginary  part  of  the  equation.   With 
a  little  algebraic  manipulation,  we  can  derive 

z  -2z„+z_  z_-z„ 

Im  ^ r: =  ll;  it 

72      h     ^v^uv 
h 

Here  we  have  a  difference  approximation  to  the  second 


~2 2 

derivative  3  z/9v   so  in  the  limit  as  h  ->  0 ,  the  assumption 

of  convergence  compels  us  to  conclude 


-  Im  z    z   =  il;   ill 

vv   V    ^v  ^uv 

as  a  second  consequence  of  the  reflection  rule. 

Before  we  go  further,  we  should  establish  that 

these  consequences  are  consistent  with  the  solution  we 

seek.   First,  notice  that  if  we  have  a  conformal  mapping, 

then  the  normal  derivative  of  ip  becomes  84;/8n  =   ^   /\z.    \ 

and  our  constant  speed  condition  says  that 
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dn 
This  verifies  the  first  consequence  of  the  rule. 

Second,  the  constant  speed  condition  must  hold  all 
along  the  top  edge  of  the  rectangle,  so  we  can  differentiate 
with  respect  to  u.   Expressing  the  condition  as 

.2 


we  have 


z  z   =  lil' 
V  V      V 


:   z   +  z  z    =    I'd)   \b 

UV  V      V  UV      ^V  UV 


If    z  (w)    is    analytic,    z      =   -iz      and    the    left-hand   side   can 
ultimate 7-Y  be  written   as    -   2    Im    (z    z      )  .      Thus   we   have 


-  Im   z    z        -   ii   ii 

V  w        ^v^uv 

which  is   just  the  second  consequence  of  the  rule.   We 

3 

might  add  that,  had  the  reflection  law  not  kept  the  h   term 

with  ^■Ji'y.'v'   ^^   would  be  imposing  the  condition 


Im  z   z   =0 
w  V 

along  T,  which  would  be  wrong. 

Although  we  have  shown  that  these  two  consequences 

are  consistent  with  our  solution,  they  do  not  appear  to  be 

as  strong  as  the  ones  generated  by  application  of  the  exact 

reflection  laws;  namely,  z  maps  bovmdary  to  boundary  and 

x^x^+  y^y   =  0 .   In  fact  niomerical  experiments  on  a  very 

simple  plasma-containment  problem  indicated  that  the 

application  of  the  inexoct  rules  does  not  determine  a 

unique  solution  to  the  difference  equations,  as  does  the 
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exact  formulation  (with  four  point  normalization) . 

Also  some  of  the  solutions  calculated  with  the  inexact 

rule  allowed  the  derivatives  z   and  z   to  become 

u       V 

non-orthoaonal ,  so  that  x  x  +y  y  was  not  zero. 
'  u  V  -^  u-'  V 

Further  experimentation  led  to  the  suspicion  that 
the  latter  fault  was  the  crucial  defect  of  the  procedure. 
The  difficulty  is  resolved  by  modifying  the  reflection  rule 
in  such  a  way  as  to  enforce  the  orthogonality  of  these 
derivati\'es .   One  can  accomplish  this  by  adding  an 
appropriate  function  of  the  derivatives .   The  new  term 
should  have  the  following  properties: 

(1)  it  should  cause  z   to  be  shifted  in  a  direction 

so  as  to  make  z   and  z  more  nearly  perpendicular, 

and  the  amount  of  the  shift  should  be  greater 

when  the  non-orthogonality,  as  measured  by 

XX  +  y  y  ,  is  greater, 
u  V   -^  u-'  V '     ^ 

(2)  it  should  be  zero  when  x  x  +  y  y   is  zero,  so 

u  V  •^u-'v        ' 

that  under  this  condition  the  term  has  no  effect 
and  the  two  consequences  of  the  reflection  rule 
will  hold. 

(3)  it  must  not  cause  instability,  or   else  the 
iterates  will  still  diverge. 

The  first  two  considerations  immediately  suggest  that 
the  term  be  proportional  to  x  x  +  y  y  ,  which,  of  course, 
depends  on  the  cosine  of  the  angle  between  z   and  s  . 
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To  determine  the  direction  in  which  z   should  be 
shifted,  we  refer  to  Fig.  6,  where  z  ,  z„ ,  and  z   denote 
the  images  of  three  consecutive  mesh  points  on  the  boundary, 
z   and  z   denote  interior  and  reflected  points  as  usual, 
Then   z   and  z   correspond  to  — ct—  and  — ^r- —  respectively. 


When  the  angle  between  z   and  z   is  less  than  tt/2  ,  as 
^  u       V 

shown,  z   should  be  shifted  in  a  direction  90°  counter^ 

clockwise  from  the  direction  of  z  -z_,  that  is,  in  the 

direction  of  i(z  -z_).   When  the  enclosed  angle  is  greater 

than  7t/2,  one  shifts  z   in  the  opposite  direction;  the 

sign  change  is  exactly  correlated  with  that  of  x  x  +  y  y  . 
^^  ^  uvu-^v 

Ultimately  one  sees  that  the  simplest  way  of 
introducing  this  shift  into  the  reflection  rule  given  above 
is  through  a  term  of  the  form 


i ( X  X  +y  y  ) 
u  V  -^  U-'  V 


Here  we  have  approximated  the  direction  of  z  -z   by  that  of 
z„-z_  because  we  don't  want  z   to  appear  in  the  right-hand 
side  of  the  reflection  rule;  also  we  put  this  in  the 
denominator  to  match  the  other  terms  in  the  formula. 
Furthermore,  we  have  taken  account  of  the  complex  conjugation 
of  z   in  the  law. 

Geometrically  it  is  clear  that  the  shift  in  z 
should  be  of  higher  order  than  z  -z^;  otherwise  we  would 
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expect  that  this  orthogonality-correcting  procedure  would 
be  unstable.   In  addition,  the  other  term  in  the  reflection 
rule  is  of  order  0(h)  =  0(Iz^-Zq|),  and  if  the  shift  is 
not  of  higher  order,  it  may  become  the  dominant  term  and 
obscure  the  free-boundary  condition.   Recall  that  this 
term  also  appears  to  lowest  order  in  the  steepest  descent 
equations,  and  that  our  goal  is  to  produce  a  faster 
convergence  process.   Finally,  the  appearance  of  a  tangential 
derivative  (z  )  in  the  formulation  should  make  us  wary  of 
instability  unless  the  term  is  of  higher  order.   These 

considerations  lead  us  to  the  conviction  that  the  correction, 

3 
as  formulated,  should  be  multiplied  by  h  ,  and  the  modified 

reflection  rule  now  reads 

-h^^^+ih.^^   ^  ,  i(x  X  +y  y  ) 

-    -   .     ^v    ^v^uv  ,  -,,3   ^  u  V  ■'u-^v 
z.  =  z«  +  +  Ah 


+     0        ^_~^o  ^-~^0 


where  A  is  a  constant  to  be  chosen  experimentally. 

The  value  of  the  constant  A  which  gives  good  results 

for  all  the  problems  considered  and  all  mesh  sizes  seems  to 

be  about  50.   However,  experimentation  has  revealed  that  this 

type  of  reflection  law,  coupled  with  an  overrelaxation 

technique  to  solve  the  equations,  is  not  very  sensitive  to 

the  particular  value  of  A  used.   Evidently  the  correction 

successfully  achieves  our  goal,  i.e.,  the  enforcement  of 

orthogonality  of  z  and  z  ,  so  that  the  solution  is 
•^      -^     u      v' 


independent  of  A , 
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Garabedian's  basic  integral  formula  quoted  at  the 
beginning  Df  this  section  from  reference  [3]  also  provides 
a  reflection  rule  for  the  axisymmetric  case,  in  theory, 
but  it  is  much  more  difficult  to  apply  the  formula  in  a 
finite-difference  scheme  because  ^    is  not  harmonic  and 
the  integral  is  not  analytic.   However,  an  obvious  modifica- 
tion of  the  two-dimensional  rule  suggests  itself. 

Instead  of 

-9iP  2  _    '^v   _  , 


we  want 


'  v' 

,2 

1    3ip  2  _   ^v    1   _  , 

y  l^vl    y 


along  F;  differentiating  this  with  respect  to  u,  we  also  want 

1  3   ,^2 


-  Im  z    z  =-.,--   ( — ) 
vv  V   2  du  y 

These  will  be  consequences  of  the  following  modification  of 
our  two-dimensional  reflection  rule: 

2 

,2  ^v  ,  ih   9   ,^v.  2 
_h   --.  +  -„5—  5—  (— . 


~7    2   511  V  '       ,  i(x  X  +y  y  ) 
y  ^  .ii.Juv-'u-'v 


z,  =  z_  +  i +  Xh"  

where  we  again  write  in  a  self-annihilating  term  which 
urges  orthogonality.   Reasoning  in  an  exactly  analogous 
manner,  we  see  that  this  rule  implies  the  conditions 
stated  above. 

We  point  out  here  that  we  have  just  arrived  at  a 
reflection  rule  by  way  of  a  heuristic  reasoning  process 
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quite  different  from  the  procedure  outlined  in  Part  2 . 
We  do  not  know  the   location  of  T,  i.e.  we  do  not  know 
F{z,z),  and  we  have  not  invoked  the  analyticity  of  the 
curve,  yet  we  are  able  to  use  reflection  to  express  the 
boundary  conditions.   This  demonstrates  the  possibility 
of  applying  the  method  to  a  very  general  class  of 
problems . 

We  conclude  that  using  the  above  approximations 
to  the  reflection  rules  in  the  vena  contracta  problem 
should  insure  the  following  three  conditions: 

(1^       7-72-  =  ^  °''  -^  7-77  =  ^ 

,2 
1  8   ,  ,  X  2       ^  , ,    1  3   ,''^v. 


(2)     -Im(z  z    )    =   -^  ir—    {\p    )         or   -Im(z   z  )  =  •55-  ,r—  (— ^) 
vv  V    2  3u  ^v  w  V    2  du    2 


(3)  X  X   +  y  y   =  0  . 

u  V    -^  u-'  V 

Condition  (3)  is  essential  for  uniqueness  for  the 
difference  equations  and  for  analyticity;  as  shown  in 

Part  1,  if  it  holds  on  the  boiindary  it  holds  everywhere, 

2   2   2   2 
and  therefore  its  harmonic  conjugate  (x  +y  -x  -y  ) /2  must 

-"^      v-'v  u-'u^ 

be  a  constant,  which  will  be  zero  when  z (w)  is  analytic. 
Equation  (1)  expresses  the  constant  pressure  condition 
when  z (w)  is  analytic,  and  (2)  is  just  a  consequence  of 
(1)  and  analyticity. 

Clearly,  using  a  finite  difference  approximation 
to  the  reflection  rule  requires  more  care  than  the  exact 
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formulation.   First,  we  can  only  conclude  that  (3)  holds 
if  our  solution  is  insensitive,  to  the  value  of  A;  this 
must  be  tested.   Second,  we  have  no  direct  proof  that  z 
will  map  boundary  to  boundary ;  our  argument  for  this  is 
based  on  the   observation  that  the  solution  to  the 
difference  equations  seems  to  be  unique,  as  indicated  by 
the  fact  that  the  overrelaxed  iterative  method  of  solution 
does  converge,  while  we  would  expect  it  to  oscillate  if 
there  were  more  than  one  solution.   And  since  our  equations 
are  certainly  consistent,  our  solution  must  be  the  right 
one,  with  the  right  boundary  correspondences.   Third,  we 

have  the  problem,  inherent  in  any  reflection  scheme  using 

2   2   2   2 
four  point  normalization,  that  x  +y  -x  -y  need  not  be  zero; 
^  '       v  •'v  u  -'u 

this  will  be  taken  care  of  in  section  F. 

To  summarize,  Garabedian's  integral  formula 

,   r       _        m/2         m/2 
(x,y;m)  =  Im  ijj^  J  (z-t)    (z-g(t)) 

^1 

-  1/2 

^(-^,-^;l;  (^-t)(z-g(t)))  ^.^^^     ^^ 

-f  ^        ^  (i-t)  (z-g(t)) 

provides  a  reflection  law  along  the  free  boundary  in  both 
the  plane  and  axisymmetric  cases.  In  two  dimensions,  the 
law  can  be  written  simply  as 


z 
f 


.2 


dz 


^0 


z+  -  Zq  "  g(z_)  -  g(zQ)  = 
but  in  three  dimensions  it  is  much  more  complicated. 
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Our  finite  difference  formula  approximating  the  rule  is 
given  for  either  case  by 

where  we  use  difference  approximations  to  the  derivatives 
appearing  on  the  right-hand  side. 

E.   The  Separation  Point 

When  the  free  boundary  reflection  rule  is  applied 
along  the  top  BC  of  the  rectangle  and  the  y-axis  rule  is 
used  on  the  side  AB,  the  corner  B  is  naturally  mapped  onto 
the  intersection,  i.e.,  the  separation  point  z  =  iY-  (Fig.  5) 
Since  the  free  surface  separates  smoothly  from  the  wall, 
the  angle  between  these  two  curves  is  180°,  while  the 
corner  at  B  is,  of  course,  90°.   This  behavior  dictates 
that  we  look  more  closely  at  the  reflection  rule  here. 
However,  as  we  mentioned  in  section  B,  the  doubling  of  this 
angle  is  desirable  for  accuracy  since  ^   has  an  expansion 
in  powers  of  (z-iY.)  '  . 

A  detailed  picture  of  the  corner  B  and  its  image 
is  given  in  Figure  7,  where  the  points  of  interest  w- , . . . ,w_ 
are  indicated.   The  corner  B  is  denoted  by  w- ,  the  neighbor- 
ing mesh  points  are  w-,  and  w^  on  AB  and  BC ,  respectively, 
and  the  neighboring  interior  point  is  w-^ .   The  reflections 
of  w^  and  Wg  through  the  side  AB  are  w^  and  w_   respectively. 
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and  w_  and  w.  reflect  w,  and  w^  through  BC.   The  z-images 
of  these  points,  carrying  corresponding  subscripts,  are 
also  shown.   The  separation  point  is  z.,   z,  lies  above  Zq 
on  the  Y-axis,  z,  is  below  z„  on  the  free  boundary  T, 
and  z^  lies  inside  the  flow  region. 

The  mapping  z (w)  doubles  the  angles  at  the  corner 
Wf.,  so  points  which  lie  180°  apart  in  the  w-plane  become 
360°  apart,  or  nearly  coincident,  in  the  z-plane,   For 
example,  w,  and  its  reflection  w-  lie  on  opposite  sides  of  Wq 
on  the  v-axis,  so  z^  and  z„  lie  on  the  same  side  of  z-  on 
the  y-axis.   Similarly  the  points  z.  and  z^,  and  z-  and  z_ 
coalesce . 

The  rule  for  reflecting  through  the  y-axis 
z^  =  -  i_ 

is  exact  and  can  be  applied  accurately  at  this  corner, 
yielding  z^  and  z_  from  z-,  and  z,.   But  our  approximation 
for  the  rule  for  reflecting  through  T   becomes  useless; 
^        and  z  ,  which  appeared  in  the  numerator  and  denominator, 
respectively,  are  both  zero  here,  because  iJj  =  1  on  AB  and 
z (w)  has  a  branch  point  at  B .   A  more  accurate  rule  must  be 
derived  for  reflecting  z^  through  z^ . 

Therefore  we  go  back  to  the  integral  representation 
in  two  dimensions: 

r 


^1 


^2  -  ^0  = 


C'^  dz  . 


^0 


As  we  indicated  above,  ?  has  an  expansion  in  powers  of 
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1/2 
{z-z_)  '     .      More  explicitly, 


?  =  Cq  +  Aj^Cz-Zq)-"-/^  +  A2(z-Zq)  +  A^Cz-Zq)  ■^/^+.  .. 
where  A^  is  zero  (cf .  [7]) .   This  gives 

so  we  can  write 


z. 


Since  z_  is  the  separation  point,  z-  =  0 .   Thus 

2  3 

z,-z„  is  0(h  )  and  the  error  here  is  0(h  ) ,  while  the 

2 

left-hand  side  is  0(h  ) .   Also  we  have 

but  analyticity  requires 


on  AB,  so 


A   =  ilj   =  0 
^x    ^y 


^0  =  i  ^x 


and  the  integral  representation  yields 

^2  -  ^0  =  -  ^^'l-^O^'^x  +  °^^^^  • 

In  order  to  use  this  as  a  reflection  rule  we  must 
find  a  sufficiently  accurate  approximation  for  ^      at  Zq . 
We  expand  i|j  around  z^  to  get  4j(x^,y^)  =  ip^,  namely 
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since  ij;  =  1  on  the  y-axis ,  ip        vanishes  and  the  above 
estimates  yield 

'''3  -  '''O  =  '''x^''3~''0^  ^   °^^^^ 
giving  the  approximation 

4;   =  — +  0(h  )  . 

Using  this  in  the  approximate  equation  for  the 
reflected  point  Z-  =  ^(^j^  '   ^®  obtain 

^3-1^0  2       3 

Therefore  we  take  as  our  reflection  law  at  the  corner  B, 

,  ,'J^3-'^0,2 
Z2  =  Zq  +  (Zq-z^)  {^^-^)       . 

Let  us  examine  the  consequences  of  this  rule. 
Since  (z_-ZQ)/h  and  (ZQ-z^)/h  are  approximations  for  3z/3v 
at  w„ ,  if  we  assume  convergence  of  the  scheme  and  take 
limits  as  h  ■*•  0  we  get 

,2 

z   =  z   ip 
V     V  ^x 

From  this  we  can  conclude  either 

(1)  Im  z  =  0    and    ^p'^   =   1    , 

V  ^x 

or 

(2)  z^  =  0  . 

But  x=OonAB,  SOX  =0  and  either  conclusion  leads  to  z  =0 , 

'      V  • 
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2 
With  this  in  mind  we  divide  by  h  /2  and  take  the 


limit.   We  get  z   =  — (z   )  ili  .   Now  we  can  conclude  Re  z   =0 
^    vv      w  ^x  vv 

2 
(which  is  no  surprise)  and  either  Im  z   =0  also  or  <iJ   =1. 
^  >^  '  vv  X 

3 
In  the  former  case,  we  have  z   =0  and  we  divide  by  h  /3' 

'  vv  J  /  r 

etc.,  etc.,  until  we  eventually  get  a  non-zero  derivative 

2 
z      ,  at  which  point  we  will  finally  conclude  ii     =1. 

W  •  'V'  '^  -^  ^x 

Therefore  we  see  that  the  consequences  of  this 
reflection  rule  are 

2 

z  =  0         and        i|^  =1. 
V  ^x 

The  first  condition  is  clearly  consistent. 
Since  ^p      is  the  normal  derivative  of  ^    in  the  z-plane  at 
the  separation  point,  the  second  equation  is  the  constant 
Speed  condition,  and  thus  it  is  consistent  too. 

Summarizing,  we  have  shown  that  in  the  two-dimensional 
case  the  reflection  rule  for  obtaining  z^   =   zCw-)  from  z^ 
should  not  be  the  same  rule  as  was  used  along  BC;  instead 
it  should  be  calculated  from 

z,  =  z„  +  (z^-zj  {— -)'^     . 

2     0      0   1   ^x^-x-. 

In  three  dimensions,  this  is  changed  to 

'^B'^O  2  1 
Z2  =  Zq  +  (Zq-z,)  (^)2  1^   . 

3   0    y 
Entirely  analogous  calculations  show  that  this  law  implies 

-T  '^x  =  1  ' 

y 

which  is  the  constant  speed  condition  at  separation. 
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In  closing  this  section  we  add  a  note  about  the 
implementation  of  this  theory.   One  desired  effect  of  this 
reflection  rule  is  to  cause  the  coefficient  of  (z.-z^), 
which  approximator,  the  speed  of  the  flow  at  z    ,    to  approach 
unity  as  h  goes  to  zero.   However,  since  the  speed  is 
already  forced  to  assume  this  value  along  the  free  boundary 
by  the  other  reflection  rule,  one  might  assume  that  this 
coefficient  may  well  be  replaced  by  1  in  the  formula, 
merely  allowing  continuity  to  enforce  the  constant  speed 
condition.   The  results  of  numerical  experiments  confirm 
this  reasoning,  and  we  found  that  comparable  accuracy  could 
be  achieved  by  the  simple  rule 

Z2  =  Zq  +  Zq  -  z^ 

in  both  the  planar  and  axisymmetric  cases . 

F.   The  Line  of  Truncation 

The  values  of  x,  y  and  i^j   on  the  vertical  line  CD 

truncating  the  infinite  rectangle  in  the  w-plane  will  approach 

their  limiting  values  at  infinity  as  the  length  of  the 

rectangle,  AD,  hereafter  called  U    ,  is  increased  (see  Fig.  5) 

max 

As  an  approximation  we  set 

y  =  v  ,      X  =  x.  =  constant, 
and 


ij,  = 
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on  CD,  and  solve  the  resulting  set  of  difference  equations 

for  several  values  of  U    ,  then  extrapolate  our  data  for 

max '  '^ 

U    -*■  °°.   To  do  this  we  must  decide  how  to  select  x,  . 
max  1 

Recall  that  no  provision  has  yet  been  made  for 

2   2   2   2 
driving  the  constant  H  =  x  +y  ~x  -y   to  zero.   Since  x,  is 
^  V  ■'v  u  -'u  1 

the  only  undetermined  parameter  left  in  the  problem,  it  is 
here  that  we  enforce  this  condition.   Experimentation  on 
several  models  has  indicated  that  the  solution  to  the 
difference  equations  is  unique  for  any  fixed  x, ,  and  further 
that  H  is  a  function  of  x, .   Therefore  it  is  through  our 
choice  of  x,  that  we  make  H  =  0 . 

If  x-  is  too  large,  the  stretching  in  the  u-direction 
of  the  map  z (w)  is  greater  than  in  the  v-direction ,  so  H 
will  be  negative;  on  the  other  hand  if  x^  is  too  small  H  will 
be  positive.   We  can  use  this  knowledge  to  calculate  x,  while 
we  are  iteratively  solving  the  difference  equations.   We  could 
calculate  H  at  each  iteration,  then  shift  x,  by  an  amount 

6  X,  =  y  H 

where  y  is  a  positive  constant.   This  is  essentially  the 
scheme  in  [1] . 

The  scheme  which  we  shall  use  is  a  slight  variation 
of  this  which  is  faster  and  more  stable.   More  specifically, 
H  is  calculated  by  finite  differences  at  every  point  in  the 
w-rectangle  to  the  right  of  u  =  1,  and  an  average  is  formed. 
Although  H  should  be  constant  once  the  solution  to  the 
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difference  equations  has  been  achieved,  it  need  not  be 
constant  during  the  iterations,  and  this  averaging  makes 
the  procedure  more  stable.   Since  we  shall  exclude  from 
our  average  the  part  of  the  rectangle  containing  the  pole 
and  the  branch  point,  the  process  is  stabilized  even  more. 

We  calculate  6x,  by  the  above  equation;  however, 
we  not  only  shift  x  by  this  amount,  but  all  the  x-coordinates 
of  the  mesh  points  to  the  right  of  u  =  1  are  shifted  by  a 
proportional  amount.   This  not  only  stabilizes  the  scheme, 
but  also  speeds  up  convergence. 

Finally,  rather  than  apply  this  procedure  at  each 

iteration,  we  wait  an  appropriate  number  of  iterations  after 

shifting  x,  for  the  effect  to  settle,  so  that  the  next 

calculation  of  H  is  more  meaningful  and  realistic.   This 

also  increases  stability.   Experimentation  indicates  that 

the  optimal  number  of  iterations  performed  between  shifts 

of  x^  should  be  about  equal  to  the  number  of  mesh  points 

along  the  u-axis ,  or  in  other  words  U   /h. 
^  '  max' 

We  can  estimate  the  optimal  value  of  the  constant  y 
by  the  following  argument.   If  a  4x1  rectangle  is  mapped  onto 
a  (4a) xl  rectangle  using  our  reflection  laws,   the  mapping 
is  given  by 

X  =  au  ,     y  =  V  . 

This  gives 

2^2     2     2    T     2 

H  =  x   +y   -X   -y   =l-a 
V    -'v     u    -^u 


and 


x,  =  4a  . 
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The  mapping  is  conformal  only  when  a  =  1,  which  means  x, 
must  be  shifted  by 

6x,  =  4  -  4a  . 

Thus  y  ought  to  be  the  ratio 

•5^1  _  4-  4a  _   4 

^  -  "T"  -  : J  -  T+E 

1  -a 

giving  y  =  2  when  a  =  1. 

We  anticipate  that  when  our  w-rectangle  in  the 
vena  contracta  problem  is  4x1,  the  optimal  value  of  y 
in  the  shift  formula  should  be  aro\ind  2.   Actual 
experimentation  shov/s  that  the  best  choice  is 

y  =  2.8  . 

G.   The  Singularities  at  the  Origin 

As  we  mentioned  in  Section  B ,  the  corner  A  of  the 
rectangle  in  the  w-plane,  located  at  the  origin,  is  mapped 
to  the  flow  region  at  infinity  in  the  second  quadrant  of 
the  z-plane  (Fig.  5) .   This  requires  a  pole  in  the 
transformation  z (w) ,  so  that  near  w  =  0 

z(w)  ^  ^ 

w 

where  R  is  a  negative  number. 

This  behavior  naturally  produces  an  unbounded 
truncation  error  in  Laplace's  difference  equation,  so  it 
must  be  compensated  for  in  some  way  before  the  problem  can 
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be  handled  on  the  computer.   The  method  we  use  is  similar 
to  that  in  [1] .   In  the  triangular  region  of  the  rectangle 
determined  by  the  base  AD,  the  side  AB,  and  the  line 
running  from  vr  =  i/2  to  w  =  1/2   (Fig.  8)  ,  we  subtract  off 
the  pole  and  solve  for  the  function 

?(w)  =  z(w)  -  ^   . 
w 

This  function  is  still  analytic;  thus  it  too  can  be 
approximated  as  a  solution  to  Laplace's  difference  equation. 
Furthermore   ?;  (w)  has  no  singularity  in  the  triangular 
region.   Therefore  the  truncation  error  is  bounded  and 
a  machine  solution  is  feasible. 

Clearly,  knowing  ^ (w)  and  the  residue  R  is  equivalent 
to  knowing  z (w) .   It  remains  for  us  to  specify  the  boundary 
conditions  to  go  with  the  difference  equation  for  ? (w) ,  and 
to  calculate  R.   We  will  discuss  the  boundary  conditions  first, 

There  is  no  problem  in  specifying  the  values  of  ? (w) 
along  the  hypotenuse  of  the  triangle,  because  we  have  values 
of  z (w)  at  these  points;  we  simply  subtract  the  values  of  R/w, 
The  details  of  the  interfacing  are  quite  straightforward, 
and  we   omit  them  here  (cf.  [1]) . 

The  boundary  conditions  for  C  along  the  two   sides 
of  the  triangle  enclosing  the  corner  A  are  derived  from  the 
reflection  laws  for  z.   Since  AB  maps  into  the  y-axis, 
the  law  for  z  is 

z^  +  z_  =  z  (w_^)  +  z  (w_)  =  0 
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where 

w=-h  +  iv,     w  =h+iv 

with  h,  as  usual,  denoting  the  rr.esh  size  and  v  denoting 
the  ordinate  of  the  mesh  point. 

Substituting  r,   +  R/w  for  z  and  noting  that  R  is 
real,  we  see  that  the  reflection  rule  for  c  is  given  by 

C+  +  ?_  =  0 

i.e.,  it  is  identical  with  the  rule  for  z. 

Similar  considerations  shov7  that  the  law  for 
reflecting  r,   along  the  bottom  of  the  triangle  is 

?+  -  C_  =  0  . 

The  procedure  for  determining  ? (w)  is  now  completely 
specified  except  for  the  calculation  of  the  residue  R. 
We  turn  our  attention  to  this  problem. 

First  we  point  out  that  if  the  function  z (w)  is 
to  be  analytic,  the  residue  R  is  determined.   We  have 
completely  specified  the  flow  region  by  setting  Y^  =  1 
(cf ,  Section  A) ,  and  we  have  fixed  the  images  of  three 
boundary  points  on  the  semi-infinite  rectangle  in  the  w-plane 
(the  point  at  infinity  and  the  two  corners) ,  so  there  is 
no  more  freedom  left  in  the  choice  of  the  confoonnal  map. 

Numerical  experiments  were  performed,  solving  the 
system  of  equations  described  above  v^ith  fixed  preassigned 
values  for  R  to  see  the  effect  of  this  parameter.   It  was 
found  that  the  iterates  converged  to  solutions  wherein  the 
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free  boundary  detached  from  the  fixed  wall  at  some  non-zero 
angle,  rather  than  detaching  smoothly  as  the  theory  predicts. 
Furthermore,  in  the  neighborhood  of  the  separation  point 
the  Cauchy-Riemann  equatioTiS  were  ivolated.   However,  we 
noticed  that  this  angle  at  detachment  was  a  monotonic  function 
of  the  parameter  R,  so  that  any  desired  angle  could  be 
achieved  by  the  proper  choice  of  the  value  of  R. 

This  immediately  suggests  an  iterative  procedure 
for  calculating  the  residue;  one  starts  with  an  initial  guess 
and  then  adds  to  this  value  a  correction  depending  on  the 
angle  formed  by  the  free  surface  and  the  wall  at  detachment, 
until  this  angle  is  zero.   The  iteration  of  the  value  of  R 
can  proceed  simultaneously  with  the  iteration  of  the 
difference  equations.   There  are  many  ways  of  implementing 
this  schem.e.   We  will  report  the  details  of  the  method 
Vvhich  seems  most  accurate,  as  judged  by  its  success  in 
the  two-dimensional  case. 

We  refer  again  to  Figure  7,  where  the  mesh  points 
along  the  top  edge  of  the  rectangle  are  called  w„ ,  w^- ,  Wq, 
and  Wq  and  their  z-images  are  similarly  designated.   The 
smooth  detachment  property  at  the  separation  point  is 
demonstrated  analytically  in  the  power  series  of  v/  along 
the  top  edge,  expanded  at  w  =  w„ .   Denoting  the  abscissa 
by  u  (the  ordinate  is  i) ,  we  have 
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x(u+i)  = 


y(u+i)  = 


^0  -^ 

8x 

u  +   — ^ 
0              8u 

^0    + 

ay 

0              9u 

u2  ^  a3^ 

^+   O(u^) 

0  ''• 

^+   0(u3)     . 

Recalling  that  x-  =  0 ,  all  first  derivatives  are  zero,  and 


we  see  that 


and 


X   =  -  X   =  0  , 

uu      w 


X  =  O(u^) 


y  -  y„  =  0(u") 


at  separation.  Hence 


y-y 


=  0(u) 


0 


so  that  the  tangent  to  the  free  boundary  is  vertical  at  z. 
and  we  have  smooth  detachment. 

Clearly  the  crucial  fact  here  is  that  the  first. 
three  coefficients  in  the  expansion  of  x  are  zero.   Our 
difference  scheme  drives  x^  to  zero  through  the  reflection 
law  on  AB,  and  x    is  driven  to  zero  through  Laplace's 
difference  equation  by  equating  it  with  -x    which,  again, 
is  zero  because  of  the  reflection  law  on  AB .   Let  us  then 
introduce  a  mechanism  for  nullifying  the  coefficient  of  u 
in  the  expansion.   This  coefficient  will  control  our 
iterations  on  the  residue  R. 

We  therefore  postulate  that  the  value  of  the  residue 
shall  be  changed  by  an  amount  proportional  to  x   evaluated 
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at  separation: 


0 

The  constant   n  is  chosen  experimentally  to  achieve  optimal 
convergence,  but  once  the  system  has  been   solved,  the 
solution  should  be  insensitive  to  changes  in  n,  since  x 
must  be  zero. 

The  best  formula  for  calculating  x   in  terms  of 
differences  involving  neighboring  mesh  points  is 


9x 


=  (3xg  -  I  Xg  +  ^  Xg)/h  +  0{h^) 


This  is  the  exact  formula  for  the  derivative  of  a  cubic 

polynomial,  so  its  use  is  most  appropriate  here,  where 

we  seek  the  dependence  in  the  form 

3 
x  '^  u 

near  separation. 

We  can  summarize  our  procedure  for  calculating  the 

residue  R  as  follows.   It  was  observed  that  the  calculated 

angle  between  the  free  boiondary  and  the  fixed  wall  is  a 

monotonic  function  of  the  value  used  for  R.    Since  the 

theory  states  that  this  angle  must  be  zero,  this  property 

is  used  to  find  the  correct  value  of  R.   An  accurate  measure 

of  the  deviation  of  this  angle  from  zero  is  provided  by 

the  derivative  x   as  approximated  through  the  above 

difference  formula,  and  this  value  is  used  in  evaluating 
the  corrective  tern 
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5R  =  n  5^ 

in  an  iterative  scheme  for  calculating  R. 

This  iterative  scheme  proceeds  along  V7ith  the 
iteration  of  the  difference  equations  for  z  and  ii , 
but  rather  than  correcting  R  at  each  cycle,  we  allow 
the  effect  of  a  correction  to  settle  down  before  correcting 
again.   In  fact,  the  iteration  of  R  is  done  simultaneously 
with  the  shifting  of  x   as  described  in  Section  C. 

We  add  that  the  above  process  works  satisfactorily 
for  the  vena  contracta,  and  that  when  the  iterations  have 
converged  the  Cauchy-Riemann  equations  hold  at  w  .   The 
cubic  approximation  for  x  gives  significantly  higher 
accuracy  at  every  mesh  size.   The  optimal  value  for  r\ 
seems  to  be  about  30 . 

We  conclude  this  section  with  a  brief  discussion 
of  the  equation  for  i|^  near  the  origin.   Since 

rP    =    1 
on  AB  and 

li;  =  0 
on  AD,  this  discontinuity  will  affect  the  accuracy  of  the 
difference  equation  unless  it  is  properly  handled. 

The  treatment  here  proceeds  exactly  as  in  [1] . 
In  the  triangular  region  we  subtract  off  a  solution  ^,, 
which  incorporates  the  behavior  near  the  origin,  so  that 
the  remainder  has  continuous  boundary  values ,   In  two 
dimensions. 
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2  ,  u  . 

ill,    =   —  arccos    ( — ,r — t-t-tt' 

1     ^  TiT+vV^ 


and  in   three   dimensions, 

I  1    ,,  u  » 

''^1    -  y  U   -  — 2 — 2   1/2'  ^ 

In  tJiree  dimensions  the  truncation  error  in  the 
difference  equation  for  ^   becomes  unbounded  at  the  origin, 
and  this,  too,  is  handled  by  a  modification  inside  the 
triangular  region.   The  variables  z  and  ip  are  transformed 
to  z'  and  ip '  via  the  Kelvin  transformation 

z'  =  1/z 

^'   =  i(;/lz|  . 

Then  the  equation  for  ^'{x',y')    has  the  same  form  as 

that  for  '|'(x,y),  and  the  truncation  error  is  again  bounded, 

For  details  see  [1] . 

H.   The  Interior  Points 

?.s  we  said  before,  the  equation 

Ax  =  0 
is  approximated  by  Laplace's  difference  equation 

X.       .     =    TT     (X.   ,   .  +  X.   .   1  +  X.,,   .  +  X.   .,,) 

1,3    T   1-1, D     i,D-l     1+1,3     ifD+1 

written  at  every  interior  mesh  point.   A  similar  equation 
is  written  for 

Ay  =  0 
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and,  in  two  dimensions,  for 

Alp  =  0  . 

In  three  dimensions,  the  equation  for  the  stream 
function,  which  we  write  in  the  conformally  invariant  form 

y 

can  be  approximated  to  the  same  degree  of  accuracy  by 
(cf.  [ID 

,|,.     .    =      i-l/J    J--1/3      1,3-1^1,3-1      1+1,3^1+1,]      i,j+l^i,3+l 

'-'  a.    ,     .+a.     .    ,+a..,     .+a.     .,, 

i-lfD      i/D-1      i+l/D      1,3+1 


where 


^'^    ^k,£+yi,j 


This  equation  is  written  at  the  interior  points  in  the 
axisymmetric  case.   After  the  Kelvin  inversion  mentioned 
in  the  last  section,  the  equation  retains  the  same  form 
for  the  transformed  values.   These  finite  difference 
equations  are  handled  by  the  standard  procedure  of 
successive  overrelaxation,  which  seems  to  be  the  best 
available  vjhen  nonlinearities  are  involved. 
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I.   THE  FORTRAN  PROGRAM 

In  Appendix  II  we  list  the  Fortran  program  used 
to  calculate  the  vena  contract  as  described  above.   A  typical 
run  involving  900  iterations  on  a  grid  of  11,000  mesh  points 
required  40  minutes  of  machine  time.   However,  due  to  the 
operating  system  on  the  CDC  6600,  the  same  result  could 
have  beer  achieved  in  about  a  third  of  this  time  if  the 
unknowns  x  =  X,  y  =  Y,  and  i|j  =  P  had  been  stored  in   single 
instead  of  double  arrays.   We  also  mention  that  when  a 
similar  large  run  was  tried  on  the  IBM  360/65  at  the 
University  cf  South  Florida  it  failed  because  of  the  smaller 
word  size. 

Careful  study  of  the  code  will  enable  the  reader 
to  identify  tJie  computations  described  in  the  text,  but 
it  will  be  helpful  if  we  define  the  input  paraur.eters . 

M,  N    are  the  number  of  mesh  points  in  the  v  and  u 

directions  respectively 
NOHALF   IS  the  number  of  times  the  mesh  is  to  be  halved 
IPOLE    is  the  location,  on  the  v  axis,  of  the  interface 

defining  the  triangle  at  the  origin  (Section  G) 
ITPSIO   is  the  number  of  iterations  performed  on  x  and  y, 

per  iteration  on  ^ 
ITERMAX  is  the  maximum  number  of  iterations  allowed 
EPS     is  m,  where  m+2  is  the  number  of  dimensions 
CON     is  X   of  Section  A,  Part  4 
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RELAXZ,  RELAXSI    are  the  parameters  c  in  the  relaxation 

2 

factors  r  =  ^ — ^-^  in  the  computation  of  z  and  ii 

respectively 
CONVERG  is  the  number  with  which  the  residuals  are  compared 

to  determine  convergence 
RESCON   is  n  of  Section  A,  Part  4 
SENS     is  y  of  Section  A,  Part  4. 


-67- 


Part  4.   RESULTS 

The  program  described  above  was  run  on  New  York 
University's  CDC  6600  computer  for  many  test  cases. 
Here  we  shall  describe  the  results  with  regard  to  the 
parameters,  accuracy,  and  convergence.   We  conclude  with 
our  new  estimates  of  the  contraction  coefficient. 

Throughout  Part  4  we  shall  compare  our  results 
with  those  of  Bloch  in  reference  [1]  ,  because  the  m.easure 
of  our  success  is  the  degree  to  which  the  reflection  scheme 
improves  on  the  efficiency  of  the  method  of  steepest 
descent.   In  Section  C  on  rates  of  convergence  the 
superiority  of  the  reflection  technique  is  clearly  demonstrated 
by  computer  runs  requiring  about  one-tenth  the  number  of 
iterations  needed  for  analogous  cases  employing  steepest 
descent. 
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A.   Dependence  on  Experimental  Parameters 

There  are  five  parameters  in  our  vena  contracta 
model  which  are  chosen  experimentally  to  achieve  certain 
goals.   These  are: 

(1)  X,      the  coefficient  of  ^^^^'''^u^v  ^^   ^^^   free 

boundary  reflection  rule,  which  enforces  orthogonality, 

2   2   2   2 

(2)  y,   the  coefficient  of  x  +y  -x^-y^  in  the  equation 

for  shifting  x, ,  the  truncated  end  of  the  jet. 

2   2   2   2 
The  goal  here  is  to  nullify  ^v+yv~Xu~yu« 

(3)  n,   the  coefficient  of  —■   at  separation,  which  occurs 

in  the  equation  for  correcting  the  residue  R  in 
order  to  achieve  smooth  detachment. 

(4)  the  location  of  the  hypotenuse  of  the  triangle 
bounding  the  region  containing  the  origin  of  the 
w-plane;  in  this  region  singular  terras  are 
subtracted  off  from  z  (w)  and  ip  (w)  .   In  the  text 

we  choose  this  line  to  run  from  w  =  i/2  to  w  =  1/2, 
but  other  choices  are  peinnissible. 

(5)  U    ,   the  length  of  the  rectangle  in  the  w-plane . 

max ' 

This  rectangle  m.ust  be  long  enough  to  produce 

accurate  approximations  to  the  theoretical  results 

of  the  infinite  jet. 

It  is  clear  that  the  rate  of  convergence  of  the 
iterates  of  the  overrelaxed  system  of  equations  will  be 
affected  by  the  choice  of  A,  y,  and n  ,  but  the  final  solution 
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must  be  insensitive  to  them  if  they  are  to  be  effective 

in  achieving  their  goals;  i.e.,  the  terms  which  they 

multiply  must  go  to  zero.   To  get  fastest  convergence, 

we  begin  the  iterations  with  the  values  A  =  50,  y  =  2.8, 

n  =  40. 

Once  convergence  is  achieved,,  as  indicated  by  all 

residuals  being  less  than  some  fixed  number,  usually  10"  , 

we  double  and  triple  each  of  the  numbers  X,  y  and  n  to 

test  the  insensitivity  of  the  solution  to  them.   In  the 

cases  tested,  this  perturbation  required  only  a  few  more 

iterations  (less  than  10%  of  the  number  of  iterations 

needed  for  initial  convergence)  to  reproduce  the  convergence 

criteria,  and  the  values  of  the  contraction  coefficient  were 

changed  by  no  more  than  one  unit  in  the  fourth  decimal  place. 

The  location  of  the  interface  between  the  triangular 

subregion  at  the  origin  and  the  rest  of  the  rectangle  should 

have  very  little  effect  on  the  data  for  small  mesh  size, 

since  the  term  being  subtracted  is  an  exact  solution  of  the 

appropriate  differential  equation.   Testing  this  witli 

different  choices  of  the  interfacing  line  produced  no  change 

in  the  fourth  decimal  of  the  contraction  coefficient  nor 

in  the  rate  of  convergence  of  the  iterates . 

The  effect  of  truncating  the  rectangle  at  u  =  U 

^  ^  max 

can  be  estimated  by  performing  the  computations  for  several 
values  of  U^  ^  and  extrapolating  for  U    •*■  °°.   In  fact, 
Bloch  [1]  has  argued  that  the  extrapolated  contraction 
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coefficient,  C  ,  differs  from  the  approximate  va3.ue 
C  (U   )   according  to  the  asymptotic  formula 

C  =  C  (U    )  +  A  e"'^  "n^ax 
c    c  max' 

where  A  and  X   are  constants .   We  shall  examine  this 

phenomenon  carefully  in  Section  D,  but  first  we  turn 

our  attention  to  the  accuracy  and  convergence  of  the  method 

on  a  rectangle  of  fixed  dimensions, 

B.   Accuracy  of  the  Solutions 

In  this  section  we  shall  verify  that  as  the  mesh  size 
h  is  decreased,  the  contraction  coef f icierit,  which  is 
representative  of  the  values  of  the  solutions  on  the  mesh 
points,  approaches  a  limit  with  an  error  varying  as  the 
square  of  the  mesh  size. 

We  anticipate  this  kind  of  behavior  because 

throughout  the  program  we  use  approximate  formulas  whose 

2 
error  are  no  bigger  than  0(h).   For  example,  our  reflection 

rules  are  all  at  least  this  accurate,  and  v/e  write  centered 

differences  or  high  order  extrapolated  differences  for 

derivatives.   The  data  presented  herein  supports  our 

2 

contention  that  the  accuracy  is  0(h  ) . 

We  tabulate  the  values  of  the  contraction  coefficient 
calculated  at  the  various  mesh  sizes  in  two  dimensions,  with 
Vx  =4  1/3. 
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c 

c 

h 

.60862 
.60968 
.60996 

t 

1/12 

1/2  4 
1/48 

Computations  show  that  this  data  is  fitted  excellently 
by  the  function 

C  (h)  =  .61005  -  .21  h^  , 
c  ' 

2 
so  we   conclude  that  the  convergence  is  0(h  )  and  the 

limit  of  C   is  .61005.   Of  course,  tliis  must  still  be 
c 

corrected  for  U    -*■    °°.      For  comrDarison,  we  quote  Bloch's 
max  -       '     ^ 

value  for  U    =4  1/2  as  .61265.   Thus  both  values  differ 
max      ' 

from  the  exact  value  .611015  by  about  .001,  and  the  accuracy 

is  comparable. 

In  the  axisymmetric  case  the  following  data  were 

computed  with  U    =4  1/3. 
^  max      ^ 


^c 

h 

.58631 
.59010 
.59103 

1/12 
1/24 
1/48 

These  numbers  are  described  by  the  function 


C^(h)  =  .59133  +  .70  h' 


and  we  again  confimnti  that  tJhe   accuracy  is  0(h)  . 
uncorrected  limit  here  is  therefore  ,59133. 
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The 


C.   Rate  of  Convergence  of  the  Iterations 

It  was  expected  that  the  iterates  of  the  reflection 
schene  would  converge  much  faster  than  those  of  the  method 
of  steepest  descent.   The  computer  rvins  confirm  this,  and 
we  present  here  a  comparison  of  two  of  our  runs  v/hich 
correspond  closely  with  the  runs  quoted  in  [1] . 

In  the  two-dimensional  case  we  chose  a  rectangle 

with 

U    =4  1/3,    h  =  1/48  . 
max      '    '  ' 

It  involves  10,192  interior  points  and  516  boundary  points, 
totaling  10,70  8  mesh  points.   The  computer  V7as  run  until 
all  residuals  were  less  than  10~  ;  the  iteration  was  started 
with  the  exact  solution.   The  number  of  cycles  required  for 
the  run  was  5  88.   For  comparison,  Bloch  quotes  a  two- 
dimensional  run  starting  with  the  same  data  and  cutting  off 

—6 

at  10   ;  it  employs 

U    =  6  ,      h  =  1/48 
max     '  ' 

and  14,161  mesh  points.   The  run  requires  6,260  iterations. 

For  the  axisymmetric  case  the  reflection  program  was 

run  with  U    =4  1/3  again,  so  that  V7hen  h  =  1/48,  the 
max       I  ^         w  /       w 

number  cf  mesh  points  was  10,70  8  as  before.   However,  to 
initialize  the  iterations  the  following  procedure  was  used. 
The  exact  solution  in  two  dimensions  was  computed  to  start 
the  procediare,  with  a  mesh  of  h  =  1/6;   after  residuals  v?ere 
reduced  to  10   ,  the  mesh  was  halved  and  the  data  interpolated 
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for  the  new  mesh  points.   This  was  repeated,  but  for 
h  =  1/48  the  iterations  proceeded  until  the  residuals 
were  reduced  to  10~  .   The  nuriber  of  iterations  required 
at  each  mesh  size  is  presented  in  the  following  table: 


Niimber 

of  Iterations 

Mesh  Size 

257 

1/6 

162 

1/12 

851 

1/24 

946 

1/48 

We  mention  that  as  h  was  decreased,  the  amount  of  over- 
relaxation  used  was  increased  much  more  rapidly  than  is 
prescribed  in  reference  [4] .   Evidently  the  nonlinearity 
inherent  in  this  problem  makes  an  analysis  of  the  optimal 
relaxation  factor  extremely  delicate. 

For  comparison,  Bloch  quotes  13,36  8  iterations 
required  to  reduce  residuals  to  10~   in  a  case  with  10,633 
mesh  points . 

We  conclude  from  this  data  that  the  reflection  scheme 
requires  about  one-tenth  as  many  iterations  as  steepest 
descent  for  the  sarr.e  degree  of  convergence. 
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D,   Estimation  of  the  Contraction  Coefficients 


As  mentioned  in  Section  A,  the  calculation  for 
the  contraction  coefficient  should  be  corrected  for  the 
truncation  of  the  rectangle  by  a  formula  (cf.  [1]) 


C  =  C  (U    )  +  A  e 
c    c  max' 


max 


This  rule  should  be  valid  at  any  mesh  size,  with  A  and  X 
independent  of  h,  at  least  when  h  is  small  enough. 

We  tested  this  behavior  by  computing  C   for  different 


values  of  U 

max 

are  as  follows ; 


The  results  for  the  mesh  size  h  =  1/24 


u 

max 

C^,  m=0 

C^,  m=l 

3^ 

.60710 

.58963 

^k 

.60968 

.59010 

7 

.61067 

.59019 

B^ 

.61069 

.59019 

where,  as  usual,  m  =  0  in  the  planar  model  and  m  =  1  in  the 
axisymmetric  case.  These  observations  fit  the  foinnula  with 
the  following  values  of  A  and  X: 


m=0 

m=l 

A 
A 

.7 
1.5 

1.2 
2.2 
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The  values  for  X   differ  only  insignificantly  from 
those  of  [1],  as  we  should  expect.   The  coefficient  A  is  of 
the  opposite  sign,  however.   This  is  probably  due  to  the 
fact  that  Bloch  uses  different  boundary  conditions  for  y 
and  ^   on  the  truncating   line  CD. 

To  estimate  the  contraction  coefficients,  we  start 
with  our  most  accurate  values,  computed  for 

U    =4  1/3,   h  =  1/48   and  residuals  less  than  10~  . 
max      /     r  / 

2 
We  first  extrapolate  these  values  for  h  -»■  0  by  the  h 

formula  of  section  B,  then  extrapolate  for  U    -»-  «>  by 

'  ^  max      -^ 

the  above.   This  leads  to  the  final  estimates 

C  =  .61106  ,    m  =  0 
c 

C^  =  .59142  ,    m  =  1 

where  at  least  the  fourth  decimal  is  significant. 
For  comparison,  Bloch  [1]  quotes 

C^  =  .61100  +  .00002 

C   =  .59135  +  .00004  , 
c  —        ' 

and  our  computations  confirm  these  values . 

In  closing,  we  add  the  following  observation. 
The  coefficient  in  the  planar  vena  contracta  is  calculated 
from  the  exact  solution  to  be 

=  .611015... 
so  our  estimate  is  high  by  .000045.   Therefore  a  prudent 
guess  for  a  five  place  value  of  C   in  the  axisymmetric  case 
is  obtained  by  subtracting  this  same  correction,  yielding 

^c  =  .591375. 
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APPENDIX  1.  ^*^^ 


PROGRAM  REFLCTR  (I NPUT , OUTPUT ) 

INTEGER  PPS 
DATA 

COMMON  X(405,5?).YM05,52).P(405,52).M. 
CWY(40,40) ,72(40.40)  . 

CN, I  POLE, PPS, II. I?. ITER,ERX, IX. JX.FRY. I Y . JY . PRS! . IS  I # JST , B6S, H, 
CITPSI,SR.GABY4.SRP.PABY4,C0N. I, J. ITPSI0.8BNS 

COMPLEX  7 

read  l.m.n.nohalt, i  pole* itpsio, i terh ax. ep^. ton . rblax? , pel  axsi. 

cconverg.rescon.sens 

READ  l.L.JUMP 
CALL  SECOND(T) 
PRINT  141, T 
PRINT  70. L 
70   rORMAT  (IH  ,  HO,*  L  ♦) 
M«l./rLOAT(M-l) 
IlaN*l 

I2=M*2 
UMAX»N*H 

1  FORMAT  (7I10/7E10.4) 
CALL  INIT 

20  PRINT  2,M,N,IP0LF. ITPSIO. I TERK AX , FPS . CON. 9il  AX? , REL AXST , rONVERa, 
CRESCON.UMAX.SENS.NOHALF 

2  FORMAT  (#0  M  N  IPOLi  ITPilO  ITERMAX  ♦, 
C  ♦  EPS  ♦/IH  ,6I10/*0  CCNTRCL  RELAX  Z  RELAX  PS 
CI*,*       CONVFRQF        RESCCN     */lH  .^Fl9.4/ 

C*0     UMAX  SENS  NOMA  F  */ 

CIH  .?E15.4.  HO) 

IWAVE=40 

ITER-0 

IPRNT»50D   $JPRNT«1    JERRMAX«100.    «SHI FT. . 11 111111 

TEST  X    in. 

ITPS1«0 

CALL  PRNTR(NOHALF, SHIFT) 

SRs2./(l.*M*RELAX7) 

QABY4a.29*SR 

SR  ■  l.-SR 

PABY432./(1 .*H*RELAXSI  ) 

SRPbI .-PABY4 

IF  (EPS.PO.O)  PABY4«.25*PaBY4 
4  ITER  »ITFR*1 

CALL  REFLEK 

CALL  INTRR 

I3»N-i 

DO  80  J«3,M 

Y( 11, J)«Y( 13, J) 

80  P( II, J)aP( 13. J) 

IF  ( ( ITER-IPRNT) .LT.O)  GO  TO  7 

IF  (TEST.LT.ERRMAX)  GO  TO  60 
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CALL  PRNTR(NOHALr, SHIFT) 

PRINT  61 
«1   FORMAT  (*    ERROR  HAS  GROWN. *} 

CALL  EXIT 
«0   ERRMAX  =  TEST 

IPRNT«IPRNT*500 

CALL  PRNTR(0. SHIFT) 
9   IF  ( ITER.LT. ITERMAX)  QO  TO  99 

CALL  PRNTR(0, SHIFT) 

CALL  EXIT 
99   TEST«AMAX1(ERX.ERY,ERSI) 

IF  (TEST.LT.CONVERQ)  GO  TO  11 

IF  (  ITER.LT. IWAVE)  GO  TO  4 

CALL  COMPMOD(SHIFT.RESCON, I  WAVE, L) 

QO  TO  4 
7   IF  ({ ITER-JPRNT) .LT.O)  GO  TO  9 

JPRNT  «  JPRNT  ♦JUMP 

CALL  PRNTRd. SHIFT) 

00  TO  9 
11   IF  (ABS(SHIFT) .LT.CONVERQ)  GO  TC  40 

CALL  COMPMOD{SHIFT.RESC0N, 1WAVE,L) 

QO  TO  4 
40   CALL  PRNTR(NOHALF. SHIFT) 

PRINT  31,  SHIFT 
SI   FORMAT  (#0   LAST  SHIFT  WAS    ♦,E10.3) 

CALL  SECOND(T) 

PRINT  141, T 

IF  {N0HALF)12,12  .13 

13  CALL  HALFER 
NOHALFsNOHALF-l 
PRINT  14 

14  FORMAT  (*0  MESH  SIZE  HAS  BEEN  HALVED.*) 
CALL  SECOND(T) 

PRINT  141, T 
141  FORMAT  (♦    TIME  IS     ♦,F10.3) 

READ  50  ,  CONVERQ.  ITERMAX, RELAXZ, SENS. RES'^QN 
50   FORMAT  ( ElO . 4 , 110 . SEIO . 4  ) 

L«L*L 

QO  TO  20 
12  CALL  EXIT 

END 


SUBROUTINE  PRNTR ( JAX, SH IFT) 
DATA 

INTEGER  EPS 
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COMMON  X(405.5?).Y(405.52).P(405.52).H. 
CWY(40,40) ,72(4n,40). 

CN, IPOLE.FPS, Il.I9.ITER,ERX. IX, JX.FRY, lY. jy.PRST, ISl « JST , PES. M, 
C1TPSI.SR.GABY4.SRP.PABY4.C0N.I,J, itpsto.ssns 

PRINT  1,  ITER.  ERX. IX, JX.ERY. lY.JY.ERSl. ISI. JST, 
CSHirT,RES 
1  FORMAT  (•0*#  I9.E90.10,  15, 15,*     X*/inX . S50 . 10 . 15. 15. ♦      V*/ 
C10X,P20.10, 15, 15.*    PSI*/ 
C10X.E20.10,10X.   ♦      SHirT*/lH  . B20 . 10 . ?flX. ♦BES* ) 

CC»(1 ./Y(2.M*1) )**(1*EPS) 

PRINT  30. CC 
SO   FORMAT  (•0    CONTRACTION  COEFFICIENT  ■   ♦.  F11.7) 

PRINT  300,  X(N*1.?) 
300  FORMAT  (•    XCEND)"   ♦,F10,6) 

RETURN 

END 


SUBROUTINE  REFLEK 
DATA 

COMMON  X(405,5?).Y{405,52).P(4  05,52).M. 
CWY(4O,40),Z2(4O.4O). 

ON, I  POLE, EPS, 1 1. I?. ITER, ERX, I  X , JX . ERY , lY. JY.PRST , TSI . JST .RES, H. 
CITPSI.SR.GABY4.SRP,PABY4,C0N, I.J. ITPilO.SBNS 

INTEOER  EPS 

COMPLEX  7.RHS 

IPMIN2=  IPOLE  -  9 

DO  1  I-IPMIN2.N 

X( I,1>hX( 1.3) 

1  Y( I,1)B-Y( 1,3) 
j»M*l 

DO  2  I"  IPMIN?,  J 
X(l,  1 )"-X{3,  1) 

2  Yd,  1)«Y(3,  I) 
IPMIN2»IPMIN2-1 
J«IPMIN2-1 

DO  3  I=J. IPMIN? 

Zb  RgS/(CMPLX(FLOAT(  I-2)»l.  )♦»-) 

X(  1,1)  »X( 1»3)  -REAL(Z) 

Y(I,1)3-Y( 1,3)  ♦A1MAQ{Z) 

7s  RES/(CMPLX(1. .FL0AT(I-2>)*H) 

X(l,  1  )  --XCS. I )  *  REAL(Z) 

3  Yd,  I)«:Y(3,  I)-AIMAO(Z) 
J«J-1 

DO  4  1=2, J 
X( I,l)aX( 1 .3) 
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V( Ijl)n-V(I,3) 
4  Yd,  I)«Y(3.I) 

IF  (EPS)  5.5.6 

5  DO  51  I«3.N 

C«-(  .5*P(  I.J-2)  -2«P(I#J-1)  ■••1.5«P(I.J))**9 

D«  (X( I*1,J)-X( I-1.J))*(X( I,J   )-X(I.J»in 
C*    (Y(  1*1.  J)-Y(  I-1#J)  )*(Y(  I,  J   >«Y(!.J-in 

Db.5*H*C0N*D 
C  ♦.2B*(( ,5*P(I*1. J-2Ul.5*P(I*l,J)-2.*P{ 1*1 ,J-1 ))*«2 

C-( .5*P( I-l, J-2>-2.*P(I-l#J-l)*l.B*P( !-l.jn**2) 

RHS  =CMPLX(C,D) 

2»  CMPLX(X(  I. J), Yd, J)) 

Z»C0NJQ(7)*RHS/(CMPLX(X(I. J-1),Y(I,J-1 ))  .7> 

X(I,J*l>cRBAL(Z) 
51  Y( I,J*1)»-AIMAQ(7) 
52   R«(P(3,M)-P(2,J)>/(X(3,M>-X(?,J)) 

ZbCMPLX(X(2.M>-X(2,J),Y(2,M)-Y(2,  J)) 
2«CMPLX(X(2,J).-Y(2,J)>  -Z*R**2 
X(2,J*1)  -    R6AL(Z) 
Y(2,J*1>»-AIMAQ(7) 
X(3,J*1)»»X(1,M) 
Y(3, J*l>»Y(l,M) 
RETURN 

6  DO  61  183, N 

D=(X{ I*l,J)-X(I-l.J)>*(X(I,J   )-X(I,J-l)) 
C*<Y<1*1,J)-Y(I-1.J))«(Y(I.J   )-Y(I,J-l)) 

PV».B*P(I,J-2)-2.*P(l, J-1)*1.5*P(I,J> 

YU»,5*(Y(I*1.J)-Y(!-1,J)) 

PUV«.25*(P( 1*1, J-2>-P{ I-l, J-2) >-(P( 1*1 , J»n-P( T-1, J-l)> 
C*.75*{P(I*1,J)-P(I-1,J)) 

YIN»1./Y(I,J) 

C«-{PV*YIN    )**9 

90  DONaCON*M 

D«.5*D0N*D-     (PIJV*PV*YU*YIN)«PV*YIN*«2 

91  RH8«CMPLX(C«D) 
2«CMPLX(X(I,J),Y(I,J)) 

2»C0NJQ(Z)  *  RHS/{CMPLX(X(I,J-1).Y(I.J»11)  -7) 
X( I, J*l)  n    REAL(7) 
61  Y(  I,J4>1)B«AIMAQ(7) 

R«(P(3.M)-P(2.J))/((X(3,M)-X(2,J))«Yep,J)> 

Z«CMPLX(X(2,M)-X(2,J),Y(2,M)-Y(2.J)) 
Z*CMPLX(X(2,J).-Y(2,J))  -Z*R**2 
X(2,J*1)  ■  REAL(7) 
Y(2, J*1>«-AIMAQ(7) 
X(3,J*1)«-X(1.M) 
Y<3,J*1>»Y(1,M) 


-89- 


PAGE 


RETURN 
END 


SUBROUTINE  INIT 
DATA 

COMMON  X(405,52).YM05.92).P(4  05.52).M. 
CWY(40,40),72(40.40), 

CN,  I  POLE, EPS, II, I?, ITER,ERX, IX, JX,eRy,IY,jy.PRST. !8IiJST,B6S.H, 
C1TPSI,SR.GABY4.SRP,PABY4,C0N. I.J. ITPSlO.ifNfi 

INTEGER  EPS 

COMPLEX  W.  ZETA,  E*  62.  8.  L,  Z 

M1«M*1 

DO  1  1-2, N 

DO  1  J«2.M1 

Us  CMPLX(FL0AT(I-?),rL0AT(J-2))*H 

IF  {I*J-4)  2,2.3 

2  X(2,2)30. 
Y(2,2)aO. 
P(2.2)>0. 

00  TO  1 

3  ZETA0l.57O79633*U 

ZETA»  .5«(CEXP(ZFTA)  -  CexP(-ZBTA)) 
ZETA*  .63661998  *CLOQ(ZETA) 
P(I,J)*  AIMA0(7ETA) 
E  ■-1.57079633*ZPTA 
E>CEXP(E) 
E2"E«E*1. 
S«-CS0RT(E2) 
G«AIMAG(S} 

IF  (Q.LT.O.)  S»CONJQ(S) 
L«CL0G(1.*S) 
OaAIMAG(L) 

IF  (Q.LT.O.)  L»C0NJG(L) 
Zs ZETA*. 63661 998* (E-S*L) 
X(I,J)«-REAL(Z) 
Y(I.J)=-A1MAG(Z)*2. 
1  CONTINUE 
N1bN*1 

X(N1.2)3X(N,2)*H 
Y(Nl,2)«n. 
DO  4  J=3.M1 
X(Nl.J)aX(Nl,2) 
Y(Nl.J)  ■  (J-2)*H 

IF  (EPS)  5,5,6 
5  PCNl.J)  ■Y(N1,J) 
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QO    TO    4 

6  P(N1,J)«       Y(N1,J)**2 
4    CONTINUE 

P(N1.2)«0. 
ID«IP0LE-4 

RES»-(l.-EPS)*.8lQ97*EPS*(-.6055) 
IF    (EPS.EQ.O)    00    TO    8 
DO    10     I«?. U 
DO    10    J«2,M1 
10       P( I,J)«.B*P( I.J) 
8    DO    7    II«3,  ID 

CALL    FLIPZ(-1. 11) 

7  CALL    FLIPSK-l.  II) 
RETURN 

BND 


SUBROUTINE    FLIPSKK.ID) 
DATA 

COMMON    X(405.52).Y(405*92)«P(405»92).H. 
CUY(40. 40). 72(40, 40)« 
CN, IPOLE.EPS, Il.I9.ITER,ERX.lX,JX,FRY.TV.Jy.PRST,TSI*JST,RES.H. 

CITPSI.SR.QABY4.SRP.PABY4.CON.I,J.ITPST0.SeNS 
INTEGER    EPS 
DO    3    Ie2.ID 
J«ID-I*2 

W»S0RT(FL0AT((l-?)**2    *(J-2)**2)) 
IF    (EPS)    2»1.2 

2  P(I,J)    «P(I,J)    ♦    .5*K*(l.-FL0AT(!-2)/W» 
00    TO    3 

1    P(I,J)sP( I,J)*    K*2*(l.-.3l830989*AC0S(«FL1ATn.2)/W)) 

3  CONTINUE 
RETURN 
END 


SUBROUTINE  FLIPZ(K. ID) 

DATA 

COMMON  XM05,5?).y(405.52),P(4  05.52).M, 

CWY(40, 40), 72(40,40), 

CN, I  POLE, EPS, U, I?, ITER,ERX, I  X, JX , FRY , I Y. JY , PRST , ISI , JST , RES, M, 

CITPSI,SR.QABY4.SRP.PABY4,C0N, I,j. ITPSTO.SCNS 

INTEGER  EPS 

COMPLEX  7,W 
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DO  1  1=2. ID 

J«1D-I*2 

W«CMPLX( ( I-2)*H, ( J-2)*H) 

7  =  CMPLX(X(  I#  J)  .  YM.  J)  )   ♦ 
X( I, J)=RFAL(7) 

y{ I. J)s  AIMAG(7) 

RETURN 

END 


K»RES/W 


90 


11 


100 


SUBR 

COMP 

COMM 

CUY(4 

CN.IP 

CITPS 
IWAV 
SOD 
DO  2 
DO  2 
XU»X 
XV"X 
YU»Y 
YV»Y 
SODb 
SOD» 
SODb 
MleM 
SHir 
IF  ( 
FACT 
DO  1 
DO  1 
X(I. 
J«M1 
SLOP 
Da-R 
IF  ( 
RESb 
D«D/ 
DO  3 
INTF 
DO  3 
ZsD/ 
X(I. 

Yd. 


OUTIN 
LEX  7 
ON  X( 

0,40) 
OLE,P 
I.SR, 
E=IWA 
■  0. 
0  I«L 
0  J«3 
(  1*1. 
(I.J* 
(  1*1, 
(I.J* 
SOD*( 
SOD/( 
SOD/F 
♦  1 

T»-SP 
ABS(S 
0R  =  1. 
00  I- 
00  J« 
J)=FA 


E  COMPMnD(SHIFT,RESCCN# IWAVS.LI 

405,52).Y(405,92)iP(405.5?).M. 

.72(40.40). 

PS.  1 1. 12. ITER.ERX. I  X , JX . f RY . I Y. JY . PRSf . ISI . JSf .RBS, H, 

GABY4,SRP.PABY4,C0N.I,w. ! TP§T 0 . SBNs 

VF*N 

.N 

.H 

J)-X(  1-1. J) 

1)-X( I. J-1) 

J)-Y( I-l, J) 

1)-Y( I, J-1) 

XU*YV)*(XU-YV)*(YU*XV)*(YU-XV) 

4.*H**2) 

LOAT( (M-?)*(N-L*1)) 

NS*SOD 

HIFT)  .GT.H)  SHIFT"SIGN(1.,SHIFT>»H 

♦SHIFT/  X(  N*i,2) 

L.U 

2. Ml 

CTOR*X( I.J) 


Ex-FL0AT(M-l)*(3.*X(S.  J)-1.5*X(4.  J)*,^3:SS*X(9.J)) 

ESC0N*SL0PE 

ABS(D) .GT. ( .09))  DbSIGN(1..D)*.0S 

RES-D 

M 

0  I«?.N 

CsMAX0(2. IPOLP-I-1) 

0  JsINTFC.Ml 

CMPLX{ FLOAT ( I -2 ) . FLOAT ( J-2 ) ) 

J)  =  X(  I, J)-REAL(Z) 
J)=Y( I, J)-AIMAQ(Z) 
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30   CONTINUE 
RETURN 
END 


SUBROUTINE  HALFER 

INTEGER  EPS 

COMMON  X(40  5,5?).Y(405.52).P(405.52).M, 

CWY{40, 40), 72(40,40), 

CN, I  POLE, EPS, 11,12, ITER,ERX,IX,JX,PRY. TY.JY.PRST, IS  I , JST , RES. H, 

CITPSI,SR,QABY4.SRP.PABY4,C0N. I , J , I TPSI 0 . SBNR 

II-lPOLE-4 

DO  22  IIIa3,  II 

CALL    FLiPZd,  III) 
22      CALL    FLIPSKI.III) 

M1»M*1 

DO    1    II«2,I1 

IsN*3-II 

DO    1    JJ«2,M1 

JsM*3-JJ 

X(2*I-2,2*J-2)«X(I,J) 

Y(2*I-2,2*J-2)«Yri,J) 
1      P(2*I-2,2*J-2)»P( I, J) 

P(2,2)s.5*(2-EPS) 

N2bN*N 

M2«2*M-1 

DO    2    I»2,N2.2 

X(I,3)=.75*X(I.4)*.375*X{I,2>-.125*X(I,6) 

Y(I,3)s.75*Y(I.4)*.375*Y(I,2)-.125*YM,6) 
2        P( I,3)s.75*P( I.4)*.375*P(I,2>-.125*P(I.6> 

DO    21    I«2,N2,2 

DO    21    J«5,M2,2 

X(  I,  J)s.75*X(  I,J-1)*.375*X(  I,u*l>-.J99*X(  T*.l-3) 
Y(  I,  J)  8,75*  Y(  I,  J-1)*,375*Y(  I .  v>l )- .  IS^^Y  t  t  ..1.3) 
21      P{  I,  J)s.75*P(  I.  J«1)*.375*P(  I,  J*1)-.125*P(  T..1-3) 
P(2,2)bO. 

M»M2 

N«N2-1 

M1»M*1 

DO    3    Ib3.N,2 
DO    3    Js2,Ml 

X(I, J)=.5*(X( I-1.J)*X( 1*1, J)) 
Y(I,J)=.5*(Y(I-1.J)*Y(I*1,J)) 
3      P(I,J)=.5*(P(I-1.J)*P( 1*1, J)) 
H8l./FL0AT(M-1) 
X(2,3)bO. 
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paQE 


S3 


789 


987 


y( 

Y( 
X( 
X( 

y( 

Y( 

P( 
p( 
p( 
p( 
IP 
II 

DO 
CA 
CA 
U 
12 
II 
DO 
IN 
DO 
X( 
Y( 
CO 

II 

DO 
IN 
DO 

y{ 

V( 

CO 
RE 

EN 


2,3) 

3,?) 

3,?) 

2.7) 

2,7) 

3,3) 

3,3) 

2,9) 

2,3) 

3,3) 

3.9) 

OLE 

=  IPO 

33 
LL  F 
LL  F 
sN*l 
=  M*2 
0«IP 

789 
TFCs 

789 
I,  J) 
I,  J) 
NTIN 
ObIP 

987 
TFCs 

987 
I.J) 
I.J) 
NTIN 
TURN 
D 


=-PES/M 

=  0. 

=-Y(2,3) 

=  0. 

=  0. 

=  .3333*(X(4,3)*X(3.2)*X(3.4)  ) 

=.95*(Y(?,3)*Y(4.3)*Y(3.2)*Y{3.4)) 

=  0. 

=.5*(2-EPS) 

=.5*(P(?.3)*P(4,3) ) 

=  0. 
=2*IP0LE-6 

LE-4 
I  I  I«3,  II 
LIPZ(-1.  Ill) 
LIPSK-I.  Ill) 


OLE 

I- 
IPO 

J- 
S.5 
=  .5 
UE 
OLE 

I- 
IPO 

J« 
=  .9 
=  .5 
UE 


2,  110.2 

LE-I-3 

3, INTFC.9 

♦  {X( I, J-1 )*X{  I, J*l)  ) 

♦  (Y( I. J-1  )  ♦YC I. J*l) ) 


3.  110.2 

LE-I-3 

2. INTFC 

♦(X( I-l. J>*X( I*l, J) ) 

♦(Y( I-l. J)*Y( I+l, J) ) 


SUBROUTINE  INTRR 

DATA 

COMMON  X(405,59).Y(405,52).P(405,52).M. 

CWY(40, 40), 72(40,40), 

CN, IPOLE.EPS. U. I?. ITER,ERX. IX, JX,ERY. lY, Jy.PRST, IS  I . JST . RES. H, 

CITPSI,SR.GABY4,SRP.PABY4.C0N. I.J. ITPSTO.SBNR 

INTEGER  EPS 
COMPLEX  W.Z 

MleM+l 

ERXsO.     SERYnO.      SL0C»1 
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N2"N*N 

0s.5*(Y(2,M)*Y(3.Ml) ) 
eRY«ABS(D-Y(2.Ml)) 
y(2,Ml)aD 
X(2.M1)«0. 

I«2 

INTrC5lP0LE-2 

DO  111  JuINTrc.M 
00  TO  1000 
ini   CONTINUE 

111  CONTINUE 
L0C«2 

DO  1  I«3,N 

INTrCsMAX0(2. IPOLF-I) 
DO  1  JelNTFCfMl 
00  TO  1000 
11  CONTINUE 

1  CONTINUE 

IF  (ITPSI)  3,2.3 

2  IF  (EPS)   21.22.21 
22  ERSI  =0. 

DO  221  1=3. N 
!NTFC=MAX0(3. IPOLE-I) 
DO  221  JbINTFCM 
00  TO  2000 
2211  CONTINUE 
221  CONTINUE 
00  TO  3 
21  ERSIaO. 

DO  211  I»3.N 

INTFCeMAX0{3. IPOLF-I) 

DO  211  JbINTFC.M 

WHYiiY(  l.J) 

A1«1./(Y( I*l,J)fWHY) 

A2sl./(Y(I-l,J)fWHY) 

A3al./(Y( I. J*1)*WHY) 

A4al./(Y( I. J-1>*WHY) 

A5ePABY4/(Al*A?*A3*A4) 

DsSRP*P( I,J)*A5*U1*P( 1*1. J)*A2*P( I-l. J)*43*P( T , J*l )*A4*P( T , J-1 ) ) 

EsABS(D>P( I.J) ) 

IF  (E.LT.ERSI )  on  TO  3001 

ERSIiE 
ISIal 
JSIaJ 
3001  P( I. J)sD 
211  CONTINUE 

3  CALL  FLIPZ(-1.  IPOLE-2) 
CALL  FLIPZ(-1.  IPOLE-3) 
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14 


(S2 


6?11 
A21 


61 


AlO 


All 


A12 


LOCaS 

MARKa 

DO  4 

INTFC 

DO  4 

IF  (  I 

CONTI 

CONTI 

IF  (  I 

CALL 

CALL 

ITPSI 

RETUR 

CALL 

CALL 

ITPSI 

IF  (F 

DO  6? 

INTFC 

DO  62 

QO  TO 

CONTI 

CONTI 

00  TO 

MARK 

DO  61 

INTFC 

DO  61 
IF  (  1 
WaCMP 
ZsCMP 
Z2(  I. 
WY(  I. 
CONTI 
MARK 
DO  61 
INTFC 
DO  61 
P(  I,  J 
MARK 
DO  61 
INTFC 
DO  61 
WHYsW 

Al  =  l. 
A2  =  l. 
A3B1. 
A4  =  l. 


IPOLE-3 
1=2. MARK 

=2*MARK-I 
J=2. INTFC 
*-J-4)  4,4.1000 

UE 

UE 

PSI)  5,6.5 

LiPZd,  IPni  E-2) 

LiPZd.  IPnLE-3) 

ITPSI-1 


LIP 
LIP 
ITP 
S) 

la 
3*M 

Ja 
200 
UE 
UE 
7 

IP 

In 

xIP 

Ja 

J-4 

X(  ( 

x(y 

)3 

)  = 

UE 
IPO 
la 

alP 

Ja 

=  7 

IPO 
la 

MAR 
Ja 
(I. 
(WY 
(WY 
(WY 
(WY 


SI(-1.  IPOLE-2) 
SI(-1.  IPOLE-3) 
SIO 

61,69,61 
3. MARK 
ARK-I 
3. INTFC 
0 


SMARK»TP0L: 


OLE-? 
2, MARK 
OLE-I 

2.  INTFC 

)  610,611,610 

I-2)*H. ( J-2)*H) 

( I, J) . Y( 1, J) )  ♦  RES/W 

1./CaBS(7) 

-AIMAQ(7)*Z2( I.J)**? 

LE-3 

3. MARK 

OLE-I 

3, INTFC 

2( I, J)*P( 1, J) 

LE-4 
3, MARK 

K-I*3 

3.  INTFC 
J) 

(  1*1. J)*WHY) 
(  I-l, J)*WHY) 
(  I, J*1)*WHY) 
( I, J-1 )*WHY) 
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A5=PABY4/(A1*A9*AS*A4) 

DaSRP*P( I, J)*A5*(A1*P( I +1 . J ) +A2*P (I -1 

E«ABS(D-P(  I> J) ) 

TO  6121 
SJSIaJ 


J)*A3*P( T,J*1)*A4*P( I. J-l) ) 


ERSI)  Gn 

SlSlBl 


IF  <e.LT 
ERSI-E 
6121  P( I,J)=D 

613  CONTINUE 
MARKbIPOLE-3 
DO  614  I«3.MARK 
INTFCalPOLE-I 
DO  614  J«3. INTFC 

614  P( I, J)=P( I,J)/72( I.J) 

7  CALL  FLIPZd*  IPOLF-2) 

CALL  FLiPZd.  IPOLF-3) 

CALL  FLIPSKI,  IPni  E-2) 

CALL  FLIPSKl.  IPnLE-3) 

RETURN 
inOO  DsSR*X(I,J)*  GABY4*(X(I*1, J) 

E«ABS{D-X( 1/ J) ) 

IF  (g.LT.ERX)  SO  TO  1001 

ERX»E 

IXbI 

JX  =  J 
inoi  X( I, J)=D 

DsSR*Y(  I,  J)-»-GABY4*(Y(  1*1.  J) 

EsABS(D-Y( I, J) ) 

IF  (E.LT.ERY)  QO  TO  1002 

JYsJ 

ERYbE 
lYal 
1002  Y(  I,J)«D 

00  TO  (ini.ll.l4>.  LOC 

2000  DsSRP*P( I, J)*PABY4*(P( I+l, J)*F(I-1 
E«ABS(D-P( I, J) ) 

IF  (E.LT.ERSI )  GO  TO  2001 

ISlBl 

ERSIbE 
JSlBj 

2001  P( I, J)5D 

QO  TO  (2001,2211.6211).  LOC 
END 


♦  X(I-l.Ji  ♦XM.J^I)  ♦X(T.,J-1)) 


+Y(I-1 . J)*V{ I . J*1)*Y< I. J-l> ) 


J)*P( I . J*1)*P( I. J-H ) 
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APPENDIX  II.   DIMENSIONAL  PERTURBATION 

In  this  appendix  we  report  a  correction  to  a 
calculation  described  in  reference  [3], 

Garabedian  describes  therein  an  entirely  different 
method  for  estimating  certain  parameters  in  axially  symmetric 
cavitational  flows  and  jets.   In  particular,  the  vena 
contracta  in  three  dimensions  is  treated  as  a  perturbation 
of  the  two-dimensional  model  with  m  as  the  perturbing 
parameter,  where  m+2  denotes  number  of  spacial  dimensions, 
as  usual.   Hence  the  technique  is  know  as  dimensional 
perturbation . 

The  results  lead  to  the  conjecture  that  the  ratio 
of  the  radius  of  the  jet  at  infinity  to  the  radius  of  the 
aperture,  Y^^yY^,  may  be  expressed  in  m+2  dimensions  as  a 
power  series  expanded  around  m  =  0 .   Thus 

Y„(m)  _  Y„(0)  _  2 

Y^^  -  Y^  -  a^m  +  a2m   +  .  .  .  . 

Regarding  m  as  a  parameter  in  the  equations  for  the  vena 
contracta,  Garabedian  is  able  to  calculate 


and,  of  course. 


Yco(-l) 

-  n 

Yo(-i) 

Y^(") 

-   1 

Y^M 

Y.(0) 

ir 

Y^m 

Tr+2 

-9  8- 


Furthermore,  the  derivative  8Y_/9m,  evaluated 
at  in  =  0  for  constant  Y^  =  1,  is  expressed  as  follows: 

1     °° 

^^0             2      \       \       ,         3+a^      ^   tan'-'-a,.    „-l           aya 
—  =   -  —  [ — ^  +   ]  tan        ^-! — 

"^      •'      •'  (1+a   )  ^  47ra  +7Tb"+Yb 


b    0 

where 

1/2 


■2 db  da 


Y 


=  2  +  (1-a  )  e   ^    -  [4a  +(l-a  )   e   '    ] 


o  ,  r  ;i  2,  ,,   2.2  -7rb/a,  1/2 
2a+[4a  +(l-a  )  e      ] 


+  2  a  log  ^ 

(l+a)"^ 

This  integral  was  evaluated  numerically  in  [3]  and  the  value 
is  quoted  as 

^  =  -  .650544  . 


3m 
Interpolating  the  above  data  with  a  cubic  polynomial  in 

5  =  m/(m+2) ,  leads  to  the  expression 

Y 

^  =  .6110  +  .4857  6  -  .1110  6^  +  .0143  6^ 
^0 

as  an  approximation  to  the  power  series.   For  three  dimensions, 

6  is  1/3  and  the  result  is 

Y^TTT-  -"^^^^ 

leading  to  the  estimate  of  the  contraction  coefficient 

,  Y»(l)   2 
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An  indication  of  the  accuracy  of  this  estimate 
is  provided  by  the  decrease  in  the  magnitudes  of  the  terms 
of  the  polynomial.   The  terms  are,  for  m  =  1, 

.4857  6  =  .1619 
.1110  6^  =  .0123 
.0143  8^   =       .0005 

and  Garabedian  arrives  at  an  error  estimate  of  1/10  per  cent. 

Using  our  larger  CDC  6600  computer  we  re-evaluated 
the  integral  for  3Y-/3m  and  found  that  a  more  accurate  value 
is 

g  Y 

x-^  =  -  .7072  +  .0003  , 
dm  — 

leading  to  the  cubic  polynomial 

^  =  .6110  +  .5279  6  -  .1110  6^  -  .0279  6^ 
0 

which  fits  the  data  at  m  =  1,  0,  and  °°  as  before. 

The  revised  computation  of  C   for  m  =  1  is  then 

Yoo(l)  2         2 
C^  =  [v  /l\]   =  .7737^  =  .5985  +  .0002  . 
c     Yq(1) 

This,  of  course,  is  larger  than  the  value  in  [3] 

and  slightly  closer  to  our  evaluation,  although  it  still 

differs  from  the  latter  by  .007.   A  glance  at  the  magnitudes 

of  the  terns  in  the  polynomial  now  reveals  that  the  error 

estimate  m.ust  also  be  revised  upwards,  for 
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.5279  6  =  .1760 
.1110  6^  =  .0123 
.0279  6^  =   .0010  . 


Consideration  of  this  data  in  the  light  of  other 
computations  of  C   leads  us  to  conclude  that  the  finite 
difference,  calculations  are  far  Hiore  reliable. 
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